SYSTEMS 

RESEARCH 

CENTER 


Supported  by  the 
National  Science  Foundation 
Engineering  Research  Center 
Program  (NSFD  CD  8803012), 

Industry  and  the  University 

Performance  Evaluation  and  Optimization  of 
Parallel  Systems  with  Synchronization 

by  L.  Gun 

Advisor:  A.  Makowski 


Ph.D.  89-3 
Formerly  TR  89-64 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

1989 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-1989  to  00-00-1989 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Performance  Evaluation  and  Optimization  of  Parallel  Systems  with 

5b.  GRANT  NUMBER 

k3jutlll  Ulll^iUlUll 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Maryland, The  Graduate  School, 2123  Lee  Building, College 
Park, MD, 20742 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

117 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


PERFORMANCE  EVALUATION  AND  OPTIMIZATION 
OF  PARALLEL  SYSTEMS  WITH  SYNCHRONIZATION 


by 


Levent  Gun 


Dissertation  submitted  to  the  Faculty  of  the  Graduate  School 
of  the  University  of  Maryland  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 
Doctor  of  Philosophy 
1989 


Advisory  Committee: 

Associate  Professor  Armand  Makowski,  Chairman/ Advisor 

Associate  Professor  Prakash  Narayan 

Associate  Professor  Lawrance  Bodin 

Assistant  Professor  Adrianos  Papamarcou 

Doctor  Alain  Jean-Marie 
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of  parallel  systems  with  synchronization 

Levent  Gun,  Doctor  of  Philosophy,  1989 

Dissertation  directed  by:  Armand  M.  Makowski 

Associate  Professor 

Electrical  Engineering  Department 


This  thesis  considers  synchronization  issues  such  as  resequencing  and  fork/join 
in  parallel  architectures.  The  discussion  is  carried  out  in  the  context  of  K  parallel 
single  server  queues  with  general  servers  where  jobs  are  subject  to  resequencing. 
Both  performance  evaluation  and  optimal  routing  problems  are  addressed  for  such 
systems. 

In  the  first  part,  Poisson  arrivals  are  assumed  to  be  randomly  allocated  to 
different  queues  according  to  a  Bernoulli  switch.  The  distributions  of  the  various 
delays  in  the  system  are  obtained  by  sample  path  arguments.  The  problem  of 
choosing  the  switching  probabilities  that  minimize  the  average  end-to-end  delay 
is  considered.  In  addition  to  obtaining  exact  results  in  some  cases,  simple  but 
accurate  approximations  are  provided  when  the  service  time  distributions  are  ex¬ 
ponential.  The  simple  form  of  these  approximations  is  then  utilized  to  solve  the 
optimization  problem  in  the  case  when  the  service  parameters  are  unknown,  and 
a  simple  stochastic  approximation  algorithm  is  proposed.  When  the  servers  are  all 
identical,  several  useful  asymptotic  results  are  obtained  as  K  increases  to  infinity. 
Various  stochastic  monotonicity  and  convexity  results  are  also  provided  for  this 
parallel  system. 


In  the  second  part,  the  dynamic  optimization  of  the  same  model  is  investigated 
under  more  general  assumptions  for  the  arrival  process.  The  resequencing  problem 
is  combined  with  a  fork/join  problem,  where  the  incoming  packets  are  broken  into 
smaller  subpackets  for  processing  at  different  queues.  The  problem  of  finding  the 
optimal  allocation  policy  that  minimizes  the  average  discounted  and  the  long-run 
average  costs  is  formulated  as  a  Markov  Decision  problem,  where  the  cost-per- 
stage  is  taken  as  the  end-to-end  delay  of  each  packet.  In  both  cases,  the  optimal 
policy  is  identified  as  the  one  that  drives  the  workload  in  each  queue  to  a  balanced 
configuration  as  quickly  as  feasible. 
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NOMENCLATURE 


a. s.  Almost  surely. 

DS  Disordering  System. 

FCFS  First  Come  First  Served. 

Geo  Geometric  distribution. 

i.i.d.  Independent  and  identically  distributed. 

IN  The  set  of  non- negative  integers. 

P(*)  Steady  state  probability  that  a  job  has  zero  resequencing  time. 
Ve  Projection  operator  onto  a  set  E. 

IR  The  set  of  real  numbers. 

IR4.  The  set  of  non- negative  real  numbers. 

RV  Random  variable. 

[x-]+  max{0,  a;},  x  E  IR. 

x  1  —  x,  0  <  x  <1. 

V  {p  E  [0, 1]*  :  Y,k=i  Pk  =  1>  ^ Pk  <  Pk,  1  <k<K} 

U  {ae[0,l]K:  E^iu*  =  1} 

U  {uE  m.K  :  ZkLi  wjfc  =  1} 

Sjj  Kronecker  symbol  defined  by  8,j  = 


r  1  i  =  j 
1 0  i^j. 


CHAPTER  I 


INTRODUCTION  AND  SUMMARY 

Recent  advances  in  computer  and  communications  technology  have  led  to 
the  proliferation  of  complex  parallel  and  distributed  system  architectures.  On  one 
hand,  these  new  systems  offer  many  advantages  over  the  conventional  systems  such 
as  resource  sharing,  reliability  and  fault  tolerance.  On  the  other  hand,  the  parallel 
and  distributed  nature  of  these  systems  pose  fundamental  problems  related  to 
the  interactions  between  different  parts  of  the  system.  For  instance,  the  ordering 
of  the  packets  is  of  great  importance  in  almost  every  application  where  data  is 
transmitted  over  several  communication  links.  This  is  also  the  case  in  distributed 
databases  or  in  computing  systems  where  computations  have  to  be  executed  in 
a  prespecified  sequence.  In  a  parallel  system  different  portions  (tasks)  of  a  job 
can  follow  different  paths  and  overtake  one  another  due  to  the  random  nature  of 
delays  in  various  parts  of  the  system.  Therefore,  some  of  these  tasks  have  to  wait 
at  the  destination  for  those  tasks  that  have  entered  the  system  earlier.  This  type 
of  synchronization  delay  is  called  the  resequencing  delay  and  may  crucially  affect 
the  performance  of  the  system.  It  is  therefore  essential  to  understand  the  effects 
of  various  system  parameters  on  this  synchronization  primitive. 

The  main  concern  of  this  thesis  is  to  provide  an  analytical  basis  for  a  better 
understanding  of  the  resequencing  constraint,  i.e.,  the  requirement  that  jobs  have 
to  leave  the  system  in  their  order  of  arrival  (see  (1.1.1)).  This  synchronization 
constraint  is  a  basic  problem  of  interest  in  many  parallel  and  distributed  computer, 
communication  and  manufacturing  systems.  Its  analysis  cannot  be  handled  by  the 
theory  of  product-form  queueing  networks,  and  its  study  therefore  provides  new 
theoretical  and  experimental  challenges. 

1.1.  MOTIVATION 

In  order  to  provide  concrete  examples,  applications  from  areas  as  diverse  as 
packet  switching,  distributed  databases  and  parallel  processing  are  briefly  con¬ 
sidered  next.  The  reader  is  referred  to  the  recent  survey  paper  of  Baccelli  and 
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Makowski  [BaM]  and  to  the  M.S.  Thesis  of  Varma  [Var.a]  for  additional  exam¬ 
ples  where  resequencing  and  other  forms  of  synchronization  constraints  naturally 
occur. 

(i)  In  packet  switching  networks  such  as  IBM’s  System  Network  Architecture 
(SNA)  [GMN]  and  the  French  PTT’s  public  network  TRANSPAC  [Dan],  mes¬ 
sages  are  formatted  into  packets  and  are  transmitted  over  an  interconnected 
network.  In  order  to  increase  the  network  utilization,  packets  belonging  to 
the  same  message  are  routed  over  different  paths  to  reach  their  destination, 
and  these  packets  may  arrive  out  of  sequence  at  the  destination  node.  Rea¬ 
sons  for  this  include  (i)  the  random  nature  of  the  delays  in  different  parts 
of  the  network  due  to  different  link  speeds,  random  path  lengths  and  vary¬ 
ing  message  sizes,  and  (ii)  the  retransmission  of  erroneous  packets  such  as 
in  Selective- Repeat  Automatic  Repeat  Request  protocols  [Sch].  In  order  to 
achieve  message  integrity,  many  communication  protocols  such  as  SNA  re¬ 
quire  First-In  First-Out  delivery.  Consequently,  the  out  of  sequence  packets 
are  stored  at  the  receiver  and  await  the  arrival  of  the  packets  that  have  been 
transmitted  before  them. 

(ii)  In  distributed  databases,  different  storage  sites  may  contain  portions  or  com¬ 
plete  replications  of  a  piece  of  data.  For  reasons  of  reliability,  centralized 
control  is  not  allowed  in  these  systems,  and  concurrency  control  mechanisms 
have  been  developed  to  preserve  consistency  [Ell,  LeL].  Due  to  random  com¬ 
munication  delays,  update  requests  such  as  read,  write  or  delete  originated 
from  different  access  sites  may  arrive  out-of-order  to  each  storage  site.  If  the 
original  order  of  the  write  and  delete  requests  are  not  preserved  in  process¬ 
ing  these  updates,  data  in  distinct  storage  sites  will  no  longer  be  replicas  of 
each  other.  Therefore,  to  ensure  consistency,  the  update  requests  must  be 
resequenced  at  each  storage  site  before  processing. 

(iii)  In  parallel  processing  machines,  different  parts  of  a  program  are  executed 
on  different  processors.  However,  data  often  needs  to  be  exchanged  between 
the  processors,  i.e.,  after  completing  the  execution  of  a  part  of  a  program  a 
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processor  may  have  to  wait  for  the  outcome  of  some  other  process  in  order  to 
continue  executing  the  rest  of  the  program.  Consequently,  the  processor  may 
have  to  wait  idle  until  the  necessary  synchronization  is  achieved. 

In  all  three  examples,  the  underlying  generic  structure  is  the  one  given  in 
Figure  1.1,  i.e.,  a  distributed  (disordering)  system  (DS)  followed  by  a  resequencing 
buffer  (RB).  In  this  figure,  the  1R+ -valued  sequences  { An ,  n  =  0, 1, . . .},  {T„,  n  = 
0, 1, . . .},  {Dn,  n  =  0, 1, . . .}  and  {S„,  n  =  0, 1, . . .}  have  the  interpretation  that 
for  all  n  =  0, 1, . . ., 

An:  Arrival  epoch  of  the  nth  job  into  the  DS,  with  Ao  =  0; 

Tn :  Delay  of  the  nth  job  in  the  DS; 

Dn:  Departure  epoch  of  the  nth  job  from  the  RB; 

Sn :  The  system  time  (or  end-to-end  delay)  of  the  nth  job,  i.e.,  Sn  =  Dn  —  An. 


Arriving  Jobs 


Out  of  Sequence  Jobs 


Resequenced  Jobs 


Figure  1.1. 

The  Generic  Structure 

The  nth  job  arrives  into  the  DS  at  time  An  and  experiences  a  delay  of  Tn  in 
the  system.  However,  due  to  the  distributed  nature  of  the  system,  the  departure 
times  {An  +  T„,  n  -  0, 1,. . .}  of  jobs  are  not  in  the  same  order  as  their  arrival 
times  {A„,  n  =  0,1,...}.  The  nth  job,  upon  leaving  the  DS,  enters  the  RB  to 
await  all  the  jobs  that  have  entered  the  system  earlier.  Only  after  all  the  jobs 
which  have  arrived  to  the  system  before  it  leave  the  DS,  does  the  nth  job  leave  the 
RB,  i.e.,  Dn  is  defined  by  the  condition 


Dn  =  .  max  {Aj  +  Tj} 

J :  Aj  <A„ 


n  =  0, 1, 


(l.i.i) 
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Throughout  this  thesis,  condition  (1.1.1)  will  be  referred  to  as  the  resequencing 
constraint.  Under  this  constraint,  the  system  time  sequence  {Sn,  n  =  0,1,...} 
satisfies  the  recursion  [BGP] 

Sn  =  max{Tn  ,  Sn-t  -  An  +  An-i}  ,  n-  1,2,...  (1.1.2) 

with  So  =  To,  and  the  resequencing  delay  Rn  of  the  nth  job  is  thus 

Rn  =  Sn-Tn ,  n  =  0, 1, . . .  .  (1.1.3) 

1.2.  LITERATURE  SURVEY 

This  section  provides  a  brief  account  of  the  literature  on  resequencing  systems. 
The  recent  tutorial  paper  of  Baccelli  and  Makowski  [BaM]  contains  a  more  detailed 
discussion  of  some  of  the  references  given  here.  Other  forms  of  synchronization 
constraints  are  also  discussed  in  [BaM],  where  some  recent  developments  in  the 
literature  are  summarized.  In  all  the  models,  unless  otherwise  specified,  the  arrival 
process  is  Poisson. 

1.2.1.  Infinite  Server  Models 

The  earliest  paper  on  resequencing  systems  is  due  to  Kamoun,  Kleinrock  and 
Muntz  [KKM]  who  studied  the  effects  of  resequencing  in  an  M/M/ oo  queue,  i.e., 
when  the  DS  is  composed  of  an  infinite  number  of  identical  exponential  servers. 
They  obtained  the  steady  state  statistics  of  the  system  time  Sn  and  of  the  number 
of  jobs  in  the  RB,  as  well  as  the  bulk  size  distribution  of  the  output  process  from 
the  RB.  The  results  of  [KKM]  were  extended  to  the  M/GI/oo  case  by  Harrus  and 
Plateau  [HaP]  by  means  of  a  similar  analysis.  Subsequently,  Baccelli,  Gelenbe  and 
Plateau  [BGP]  considered  the  situation  where  the  RB  is  followed  by  a  single  server 
queue  with  a  general  service  time  distribution.  The  DS  is  again  implemented  by  the 
M/GI/oo  queue.  These  authors  gave  recursive  formulas  of  the  type  (1.1.2)  for  the 
end-to-end  delay  (including  the  delay  in  the  single  server  queue),  which  are  then 
used  to  derive  integral  equations  for  its  distribution.  The  analysis  is  carried  out 
via  factorization  methods  when  the  delay  distributions  in  the  DS  are  exponential, 
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i.e.,  when  DS  is  the  M/M/oo  system.  Recently,  Gelenbe  and  Stafylopatis  [GeS] 
considered  the  model  of  [HaP]  under  three  partial  ordering/resequencing  disciplines 
which  reflect  the  random  association  or  locality  in  the  precedence  constraints. 

1.2.2.  Finite  Server  Models 

In  the  infinite  server  models  of  Section  1.2.1,  the  disordering  delays  {Tn,  n  = 
0,1,...}  form  an  independent  and  identically  distributed  (i.i.d.)  sequence  of  ran¬ 
dom  variables  (  RVs)  which  is  also  independent  from  the  arrival  time  sequence 
{A„,  n  =  0,1,...}.  When  the  number  of  servers  is  finite,  these  independence 
properties  are  no  longer  available,  and  the  corresponding  results  are  few  in  the 
literature  and  limited  to  very  special  cases.  All  models  assume  Poisson  arrivals 
and  exponential  service  time  distributions.  Two  classes  of  models  are  considered 
for  the  DS,  namely  (i)  models  where  there  is  a  common  buffer  attended  by  parallel 
servers,  (ii)  models  with  parallel  single  server  queues. 

(i)  Common  Buffer  Models:  The  models  in  this  class  all  lead  to  a  Marko¬ 
vian  analysis.  However,  the  lack  of  independence  causes  a  rapid  explosion 
in  the  size  of  the  state  space  with  the  number  of  servers.  The  M/M/K 
model  with  identical  servers  is  studied  by  Bharat-Kumar  and  Kermani  [BKK] 
who  derived  an  expression  for  the  mean  resequencing  delay.  Agrawal  and 
Ramaswamy  [AgR]  also  considered  the  model  in  [BKK]  but  focused  on  the 
distributional  aspects  of  the  resequencing  delay.  Yum  and  Ngai  [YuN]  consid¬ 
ered  the  more  general  M/M/K/B  model  with  heterogeneous  servers  when  the 
common  buffer  size  B  is  finite.  In  the  case  of  heterogeneous  servers,  it  is  com¬ 
monly  assumed  in  the  literature  that  when  more  than  one  server  is  available, 
jobs  are  scheduled  to  the  fastest 'available  server.  Yum  and  Ngai  obtained  the 
resequencing  delay  distribution  through  a  numerical  algorithm.  However,  due 
to  the  high  dimensionality  of  the  state-space,  they  reported  their  method  to 
be  limited  to  K  <  5.  Lien  [Lie]  considered  the  resequencing  delay  due  to  the 
M/M/2  queue  and  obtained  the  average  resequencing  delay  by  a  simple  yet 
clever  argument.  Later  on,  Varma  [Var.a]  showed  that  the  Markovian  state 
representation  of  Lien  naturally  leads  to  a  matrix-geometric  representation  for 
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the  steady  state  buffer  occupation  probabilities  when  the  M /M/2  queue  has  a 
finite  buffer.  However,  this  approach  does  not  extend  to  the  case  K  >  2.  The 
M/M/2/B  model  is  also  considered  in  the  Ph.D.  Thesis  of  Iliadis  [Ili]  where 
the  distribution  of  the  resequencing  delay  is  obtained  under  various  threshold 
type  scheduling  policies  (see  also  [IL.a]  and  [IL.b]). 

(ii)  Parallel  Buffer  Models:  The  work  on  this  class  of  models  is  very  re¬ 
cent.  Jean-Marie  [JeM.b]  considered  the  case  of  two  M/M/1  queues  in  parallel 
with  identical  servers  when  two  different  Poisson  streams  of  incoming  jobs  are 
routed  randomly  to  the  queues  by  two  Bernoulli  switches.  He  obtained  the 
distribution  of  the  resequencing  delay  by  a  sample  path  argument.  Iliadis  and 
Lien  [IL.c]  considered  two  parallel  heterogeneous  M/M/1  queues  with  two 
types  of  traffic,  namely  (i)  jobs  allocated  to  the  queues  by  a  Bernoulli  switch 
which  are  subject  to  resequencing,  and  (ii)  interfering  local  traffic  to  each 
queue  which  are  not  subject  to  resequencing  and  leave  the  system  as  soon  as 
they  are  serviced.  All  arrival  processes  are  assumed  to  be  Poisson  and  mutu¬ 
ally  independent.  A  recursive  method  is  proposed  for  obtaining  the  average 
resequencing  time  for  the  jobs  allocated  by  the  Bernoulli  switch.  Their  solu¬ 
tion  method  extends  to  the  situation  where  there  are  several  arrival  processes 
to  the  Bernoulli  switch  from  various  sources  when  resequencing  operates  class 
by  class,  i.e.,  when  jobs  from  a  given  source  need  to  wait  only  for  the  jobs 
that  arrive  to  the  DS  from  that  source. 

1.2.3.  Structural  Results 

Bounding  methodologies  based  on  convex  and  strong  stochastic  ordering  ar¬ 
guments  are  used  in  the  M.S.  Thesis*  of  Varma  [Var.a].  He  provides  monotonicity 
results  for  DSs  in  (i)  and  (ii)  above  under  more  general  assumptions  on  the  service 
and  arrival  distributions.  Multistage  DSs  in  tandem  are  also  considered  in  [Var.a] 
where  it  is  shown  that  for  infinite  server  DSs  resequencing  the  jobs  only  after  the 
last  stage  yields  a  lower  end-to-end  delay  than  when  jobs  are  resequenced  after 
every  DS.  When  the  DSs  are  GI/GI/K  queues  and  the  jobs  are  resequenced  after 
every  stage,  Varma  also  showed  that  a  decrease  in  service  times  in  one  of  the  stages 
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implies  a  decrease  in  the  end-to-end  delay. 

1.2.4.  Optimization  of  Resequencing  Systems 

Results  on  the  optimal  design  of  resequencing  systems  are  again  very  few 
and  recent.  The  static  optimization  problem  of  choosing  the  optimum  switching 
probabilities  in  the  multiple  queue  model  of  Section  1.2.2(h)  is  considered  by  Jean- 
Marie  [JeM.b]  when  K  =  2  with  identical  exponential  servers.  Two  Poisson  input 
processes  are  randomly  allocated  to  the  queues  by  two  Bernoulli  switches  when 
the  resequencing  operates  class  by  class.  Mean  queueing  and  resequencing  times 
are  obtained,  and  the  optimal  switching  probabilities  are  characterized  through 
the  unique  solution  of  a  fifth  order  polynomial. 

Varma  [Var.b]  considered  the  dynamic  allocation  of  jobs  to  servers  in  the 
M/M/2  queue  with  heterogeneous  servers  so  as  to  minimize  the  system  time. 
Using  dynamic  programming  arguments  he  showed  that  the  faster  server  should 
always  be  kept  busy,  while  the  optimal  policy  that  assigns  jobs  to  the  slower  server 
is  of  threshold  type  in  the  number  of  jobs  in  the  common  buffer,  and  is  independent 
of  the  number  of  jobs  in  the  RB. 

In  addition  to  the  papers  annotated  above,  the  effect  of  resequencing  on  more 
general  systems  are  considered  in  [Bac]  and  [JeM.aj.  While  [JeM.a]  studies  the 
effect  of  the  resequencing  delay  in  Omega  networks,  [Bac]  considers  synchronization 
issues  in  distributed  databases  and  provides  computable  bounds  for  various  delays. 

1.3.  SUMMARY  OF  THE  THESIS 

In  this  thesis  we  consider  the  performance  evaluation  and  the  optimal  routing 
of  jobs  in  the  resequencing  system  of  the  type  given  in  Figure  1.1.  Attention  is 
given  here  to  DSs  composed  of  K  parallel  single  server  queues  with  infinite  capacity 
buffers.  The  material  in  Chapters  II  and  IV  is  joint  work  with  Dr.  A.  Jean-Marie 
and  is  also  included  in  the  manuscripts  [GJM]  and  [JMG],  respectively. 

In  Chapter  II,  jobs  arrive  according  to  a  Poisson  process  and  are  randomly 
allocated  to  the  queues  by  a  Bernoulli  switch.  Each  queue  is  attended  by  a  sin¬ 
gle  server  with  a  general  service  time  distribution.  The  service  times  at  different 
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queues  are  assumed  independent  with  possibly  different  distributions,  i.e.,  the  DS 
is  composed  of  K  parallel  M/GI/1  queues.  We  first  express  the  resequencing  and 
system  times  in  terms  of  the  waiting  and  response  times  in  each  one  of  the  parallel 
M/GI/1  queues.  This  is  done  by  means  of  a  sample  path  argument  which  focuses 
directly  on  the  system  time  rather  than  on  the  resequencing  delay.  The  distribu¬ 
tions  and  first  moments  of  the  resequencing  and  system  times  are  then  obtained 
using  these  expressions.  The  problem  of  choosing  the  switching  probabilities  in 
order  to  minimize  the  average  end-to-end  delay  is  then  addressed  for  the  following 
two  special  cases:  (i)  when  the  service  time  distributions  are  identical,  and  (ii) 
when  K  —  2  and  the  service  times  are  exponential.  In  the  first  case,  it  is  shown 
that  the  equal  load  allocation  is  optimal,  and  that  at  this  optimum  configuration 
the  average  system  time  decreases  with  K.  The  second  case  differs  from  the  model 
in  [JeM]  mentioned  in  Section  I.2.2(ii)  in  that  now  there  is  only  one  arrival  pro¬ 
cess  and  the  servers  have  different  rates.  The  optimal  routing  probability  is  again 
characterized  by  the  unique  root  of  a  fifth  order  polynomial  for  certain  values  of 
the  system  parameters,  whereas  it  is  best  to  route  all  the  traffic  to  the  faster  server 
for  other  values.  The  equations  that  define  these  regions  in  the  parameter  space 
are  identified. 

The  results  obtained  in  the  special  case  (ii)  lead  to  an  asymptotically  exact 
approximation  for  K  (>  2)  heterogeneous  exponential  servers.  This  approximation 
has  a  very  simple  form  and  provides  insights  into  the  variation  of  the  optimal 
switching  probabilities  with  the  system  parameters.  Furthermore,  for  wide  range 
of  system  parameters  it  yields  a  relative  error  of  at  most  1%  for  the  mean  system 
time. 

The  computation  of  the  optimal  routing  probabilities  for  the  same  system 
without  resequencing  is  also  recalled  and  improved.  The  solutions  to  both  prob¬ 
lems  are  compared  throughout  Chapter  II.  Strong  stochastic  convexity  (see  Defi¬ 
nition  A.II.l  of  Appendix  II)  of  the  stationary  disordering  delay  in  the  switching 
probability  vector  is  also  established.  In  the  homogeneous  case,  equal  load  alloca¬ 
tion  is  shown  to  stochastically  minimize  the  disordering  delay.  Furthermore,  at  this 
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optimum  configuration,  stochastic  monotonicity  and  integer  convexity  properties 
in  K  are  established  for  the  disordering  delay. 

Chapter  III  considers  the  model  of  Chapter  II  when  the  service  times  are  expo¬ 
nentially  distributed,  but  when  both  the  service  and  the  arrival  rates  are  unknown. 
The  simple  form  of  the  approximations  given  in  Chapter  II  is  used  in  a  stochastic 
approximation  algorithm  for  computing  the  approximate  optimal  switching  prob¬ 
ability  vector.  The  algorithm  starts  with  an  initial  switching  vector  and  estimates 
the  system  parameters  by  making  idle  time  measurements.  The  switching  vector 
is  then  updated  using  these  measurements.  The  algorithm  proceeds  by  taking  new 
measurements  when  the  system  is  driven  with  this  new  switching  vector,  and  then 
updates  the  switching  vector  again. 

Chapter  IV  studies  the  asymptotic  behavior  of  various  delays  in  the  model  of 
Chapter  II  when  the  service  time  distributions  are  identical  and  the  load  is  equally 
(optimally)  allocated  to  the  queues.  Asymptotic  expressions  for  the  distributions 
of  the  resequencing  and  system  times  are  provided  as  K  increases  to  infinity.  Two 
cases  are  considered  depending  on  whether  the  arrival  rate  into  the  system  is  held 
constant  or  grows  linearly  with  K. 

In  the  first  case,  it  is  proved  that  the  statistics  of  the  system  tend  to  those  of 
the  M/GI/oo  system  with  resequencing.  While  asymptotic  stochastic  monotonic¬ 
ity  and  convexity  results  (in  the  sense  defined  in  Chapter  IV)  are  stated  for  the 
system  time  RV,  the  resequencing  delay  exhibits  different  structural  characteris¬ 
tics  depending  on  the  arrival  rate  into  the  system.  For  example,  when  the  service 
times  are  exponential  it  is  (asymptotically)  stochastically  increasing  and  concave 
in  K  for  smaller  loads  on  the  system,  while  for  larger  loads  it  is  (asymptotically) 
stochastically  decreasing  and  convex. 

The  situation  is  different  in  the  second  case.  Indeed  the  expected  system  and 
resequencing  times  both  grow  as  log  K  while  the  average  response  time  is  constant 
for  all  K.  Therefore,  although  the  response  time  of  a  job  in  the  queue  dominates  the 
resequencing  time  in  the  first  case,  the  resequencing  delay  dominates  the  response 
time  and  has  a  major  impact  on  the  system  time  in  the  second  case. 
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In  Chapter  V,  it  is  assumed  that  the  jobs  arrive  according  to  a  renewal  se¬ 
quence.  The  workloads  of  the  jobs  are  i.i.d.  and  independent  from  the  interarrival 
times.  Upon  arrival,  instead  of  the  Bernoulli  allocation  as  in  the  previous  chap¬ 
ters,  each  job  is  broken  into  smaller  tasks  for  processing  at  different  queues.  These 
tasks  are  later  assembled  in  the  RB,  and  a  job  leaves  the  system  only  after  all  of 
its  tasks  are  serviced  and  all  the  other  jobs  that  arrived  into  the  system  before  it 
are  also  assembled  and  have  left  the  system.  The  dynamic  optimization  problem 
of  allocating  each  job  into  the  parallel  queues  is  considered  when  the  processing 
rates  in  different  queues  are  identical  and  constant,  i.e.,  the  service  time  of  a  task  is 
proportional  to  its  workload  and  tasks  with  equal  workloads  require  equal  amounts 
of  service  time  at  different  queues.  The  average  finite  horizon,  long-run  average 
and  average  discounted  costs  are  all  shown  to  be  minimized  by  the  same  allocation 
policy  when  the  cost-per-stage  is  taken  to  be  the  system  time  of  a  job.  At  each 
arrival  epoch,  given  the  workload  in  each  queue  and  the  workload  of  the  arriving 
job,  this  optimal  policy  allocates  the  workload  of  the  job  to  the  queues  so  as  to 
derive  the  workload  in  each  queue  into  a  balanced  position  as  fast  as  feasible.  A 
simple  algorithm  that  computes  this  optimal  allocation  vector  is  also  provided. 
Under  this  optimal  policy,  the  steady  state  system  time  of  a  job  is  obtained  as  the 
response  time  in  a  GI/GI/1  queue. 

The  case  where  the  jobs  are  not  allowed  to  be  broken  into  smaller  tasks  is  also 
briefly  considered  when  K  =  2.  Optimality  of  joining  the  queue  with  the  smallest 
workload  is  established  in  this  case. 

Notation 

The  following  notation  and  definitions  are  used  throughout  this  thesis:  The 
(resp.  positive)  real  line  is  denoted  by  (resp.  H+)  H.  The  kth  component  of  a 
column  vector  x  in  IRA"  is  denoted  by  a:*,  while  the  kth  component  of  the  vector 
[x]+  is  defined  as  [rr*]+  :=  max{0,£jfc},  1  <  k  <  I< .  The  transpose  of  x  is  denoted 
by  xT .  The  column  vector  of  ones  with  appropriate  dimensions  is  denoted  by 
e.  The  notation  Ve{%)  denotes  the  projection  of  x  onto  a  set  E  C  IRK.  The 
distribution  function  of  a  RV  X  will  be  denoted  by  the  same  letter,  e.g.,  X(x)  = 
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P(X  <  x).  Whenever  it  is  necessary  to  explicitly  indicate  the  dependence  of  A'(a:) 
on  a  parameter  9,  the  notation  X(9,  x)  and  AT(0)  will  be  substituted  for  X(x)  and 
X ,  respectively.  The  distinction  between  the  distribution  function  X(x)  and  the 
RV  X(9)  will  always  be  clear  from  the  context.  Finally,  3T  =  1  —  x  for  every  x  in 
[0,1]. 
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CHAPTER  II 


RESEQUENCING  IN  PARALLEL  QUEUES 
WITH  BERNOULLI  LOADING 

II.1.  INTRODUCTION 

A  system  with  K  parallel  queues  each  with  an  infinite  capacity  buffer  is  con¬ 
sidered  when  the  jobs  are  subject  to  resequencing.  Jobs  arrive  according  to  a 
Poisson  process  and  are  allocated  to  the  queues  by  a  Bernoulli  switch.  The  service 
times  at  different  queues  are  assumed  independent  with  possibly  different  distri¬ 
butions.  After  completing  service,  a  job  moves  to  the  RB  and  awaits  the  service 
completion  of  the  jobs  which  have  arrived  into  the  system  before  it.  The  model  is 
made  precise  in  Section  2. 

In  Section  3,  the  distributions  of  the  resequencing  and  system  times  are  com¬ 
puted  in  terms  of  the  distribution  functions  of  the  waiting  and  response  times  in 
the  parallel  M/G//1  queues.  Expressions  for  their  first  moments  are  then  easily 
derived. 

The  problem  of  choosing  the  allocation  probabilities  to  minimize  the  average 
system  time  is  addressed  in  Section  4  in  two  special  cases,  namely  an  arbitrary 
number  of  identical  servers  and  two  queues  with  exponential  servers.  Section 
5  contains  simple  approximations  to  the  optimal  switching  vector  when  the  K 
parallel  servers  are  exponential. 

The  solution  to  the  problem  without  the  resequencing  constraint  is  briefly 
recalled  [BuC]  in  order  to  provide  a  comparison  to  the  resequencing  problem.  An 

r-' 

algorithm  to  compute  the  optimal  probabilities  is  available  for  the  general  case  in 
[BuC].  Simplifications  to  this  algorithm  are  provided  in  Appendix  I. 
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II.2.  THE  MODEL 


The  model  consists  of  K  queues  in  parallel  where  each  queue  fc,  1  <  k  <  K, 
has  an  infinite  capacity  buffer  and  is  attended  by  a  single  server  whose  service  time 
distribution  Bk{ •)  has  finite  mean  1/fik  and  variance  The  service  times  are 
assumed  mutually  independent.  Jobs  arrive  into  the  system  according  to  a  Poisson 
process  with  parameter  A,  and  join  the  kth  queue  with  probability  pfc,  1  <k<K, 


Figure  II.  1. 
The  Model 
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where  pk  =  T  The  service  discipline  at  each  queue  is  first  come  first  served 

(FCFS),  and  upon  service  completion,  a  job  joins  the  RB  to  await  for  the  service 
completion  of  all  the  jobs  that  have  entered  the  system  before  it. 

Throughout  this  chapter,  the  system  is  considered  in  statistical  equilibrium. 
For  each  queue  k,  1  <  k  <  K ,  set 

Bk-  the  service  time  of  a  job  in  server  k ; 

Wk'  the  waiting  time  of  a  job  in  queue  k\ 

Tk‘.  the  sojourn  time  of  a  job  in  queue  k  (i.e.,  Tk  =  Wk  +  Bk), 

and  for  the  system,  define 

T:  the  response  time  of  a  job  (i.e.,  the  time  needed  to  complete  its  service); 
R:  the  resequencing  time  of  a  job  (i.e.,  the  time  it  spends  in  the  RB); 

S:  the  system  time  of  a  job  (i.e.,  S  =  R  +  T). 

Since  the  thinning  of  a  Poisson  process  by  a  Bernoulli  process  results  in  a 
set  of  independent  Poisson  processes,  the  system  behaves  like  a  collection  of  K 
independent  M/GI/ 1  systems  with  arrival  rates  Ajq , . . . ,  XpK,  and  with  a  global 
ordering/resequencing  mechanism.  In  particular,  the  RVs  Wk  and  Tjt  are  inde¬ 
pendent  of  the  RVs  Wi  and  T;  for  1  <  k  ^  l  <  K.  Since  the  arrival  process  to 
each  queue  is  Poisson,  the  RV  Wk  is  the  workload  of  queue  k  at  a  random  instant 
[Kle.a],  i.e,  Wk  is  the  virtual  waiting  time  in  queue  k. 

The  presence  of  resequencing  does  not  affect  the  stability  condition  of  the 
system  [BGP].  Therefore,  the  system  is  stable  if  and  only  if  each  queue  is  stable, 
which  reads 

0  <Pfc<y,  1  <k<K.  (2.2.1) 

The  utilization  of  queue  k  is  denoted  by  pk  =  A Pk/pk,  while  the  system  capacity  is 
p  =  Pk-  Note  that  the  stability  condition  (2.2.1)  implies  A  <  p.  Conversely, 

if  this  condition  is  satisfied,  then  the  convex  subset  V  of  IR;<"  defined  by 

K 

V  =  {p  e  [0, 1]^  :  ~  1  an<*  ^ Pk  <  Pk'  1  ^  &  <  K) 

*= l 
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is  nonempty. 


II. 3.  RESPONSE  AND  SYSTEM  TIME  ANALYSIS 

In  this  section,  expressions  for  the  RVs  T,  R  and  S  are  provided  in  terms  of 
the  RVs  Tk  and  Wk,  l  <  k  <  K.  These  expressions  then  lead  to  the  computation 
of  their  distribution  functions  and  first  moments. 

On  an  underlying  common  probability  space,  define  a  RV  U  :  Q  — >  {1,  ...,K} 
such  that  P(U  =  k)  =  pk,  1  <  k  <  K.  The  RV  U  is  assumed  independent  of  all 
the  RVs  introduced  so  far. 

Theorem  II. 3.1.  The  RVs  T,  R  and  S  are  given  respectively  by 

T  =  Tv  ,  (2.3.1a) 

R  =  \W*  -  Tv]+  (2.3.1  b) 

and 

S  =  ma x{W*,Tu}  ,  (2.3.1c) 

where 

W*  :=  max{W*;  1  <  k  <  K,  k  ^  U}  .  (2.3.1  d) 

Although  a  transient  version  of  (2.3.1)  also  holds  for  each  job  n,  n  =  0,1,..., 
only  the  steady  state  equalities  are  given  here  in  order  to  keep  the  notation  to  a 
minimum. 

Proof.  The  expression  for  T  is  plain  and  expresses  the  Bernoulli  allocation  of 
the  jobs. 

In  order  to  compute  R,  assume  that  a  “tagged  job”  C  arrives  to  the  Bernouilli 
switch  at  time  t  =  0.  Let  Ck  be  the  last  job  in  the  FCFS  queue  k  at  time  t  —  0+  so 
that  C  is  Cu •  Then,  with  the  notation  introduced  in  Section  2,  the  job  Cjt,  k  ^U, 
completes  service  at  time  Wk,  so  that  W*  is  the  time  when  the  last  of  the  CVs, 
k  7^  U,  completes  service.  The  tagged  job  C  completes  service  at  time  Tjj-  If 
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Tu  <  W* ,  then  at  least  one  of  the  Ck,  1  <  k  ^  U  <  7i",  has  not  yet  completed 
service  and  C  has  to  wait  until  W* ,  i.e.,  W*  —  Tu  is  his  resequencing  time.  On 
the  other  hand,  if  Tu  >  W* ,  then  C  experiences  no  resequencing  delay.  Therefore, 
R  =  [W*  -  Tu}+. 

The  expression  for  S  now  follows  immediately. 

□ 


The  following  lemma  presents  an  equivalent  form  of  Takacs’  integro- 
differential  equation  for  the  M/GI/1  system  [Kle.a,  p.  230]  and  will  prove  useful 
in  computing  the  distribution  of  S.  The  term  A pk  is  the  arrival  rate  into  the  kth 
queue. 

Lemma  II. 3. 2.  The  equality 


1  dWu 

T*(*)  =  Wk(x)  -  1  <k<K  (2.3.2) 

holds  for  every  x  >  0  where  Wk(x)  is  differentiable. 

The  central  result  of  this  section  is  the  following  theorem. 

Theorem  II. 3. 3.  The  distribution  functions  T(-),  R(-)  and  S(-)  are  given  by 

K 

T(x)  =  YlPkTk^  ’  (2.3.3a) 

t=i 


and 


K  foo  K 

R(x)  =  Ep*  /  +  t)dTk(t ) 

k=l  '=i 

i*k 


(2.3.3 b) 


S(x) 


K 

n  w*(x'> 


k=  1 


ld_ 
A  dx 


K 


n  wk(x) 


\k= 1 


(2.3.3c) 


for  every  x  >  0  where  Wk(x ),  1  <  k  <  K,  are  all  differentiable.  Their  means  are 
given  by 


K 


ET  = 

Jfc  =  l 


2  -  pkiff  -  a\ii%) 
2(/i k  A Pk ) 


(2.3.4a) 
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ER  =  ES-ET 


(2.3.46) 


and 


1*00  K  I  K 

ES  =  /  (!  -  II  +  T(!  -  TL1  -  «))  •  (2.3.4c) 

•'»  tei  A  tei 


Proof.  Conditioning  (2.3.1a)  on  the  value  of  U  gives  (2.3.3a),  from  which 
(2.3.4a)  follows  by  a  simple  application  of  the  Pollaczek-Khinchin  mean  value  for¬ 
mula  [Kle.a,  p.  190]. 

Equation  (2.3.36)  follows  from  (2.3.16)  by  conditioning  on  U,  and  noting  that 
P([A]+  <  x)  =  P(^4  <  x)  ,  x  >  0  , 
and  (2.3.46)  obvious  since  S  =  T  +  R. 

To  obtain  the  distribution  of  S,  note  that  the  RVs  in  the  right-hand  side  in 
(2.3.1c)  are  conditionally  mutually  independent  given  U .  The  computation,  which 
makes  use  of  Lemma  II. 3. 2,  proceeds  as  follows:  For  all  x  >  0  where  the  functions 
Wk ,  1  <  k  <  K,  are  all  differentiable,  it  is  plain  that 


K  K 

Six)  =  ^pfcTfc(x)fJkK,(^) 

fc=i 


i=i 


=  Y,pk^Wk^~ 


&=1 

K 


i=i 

l*k 


fc=i  '  it=i  <=i 


and  this  rewrites  as  (2.3.3c).  Equation  (2.3.4c)  then  follows  by  routine  integration 
since  ES  —  /0°°(1  —  S(x))dx  and  W*(0)  =  1  —  />*,  1  <  k  <  K. 

□ 

Comparison  of  the  formulas  for  T(®)  and  S(x),  or  for  ET  and  ES ,  quickly 
reveals  how  the  resequencing  requirement  complicates  the  analysis.  In  particular, 
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whereas  ET  is  expressed  in  terms  of  the  first  two  moments  of  the  service  times, 
the  distributions  of  the  waiting  times  are  needed  to  compute  ES.  This  point  is 
made  more  apparent  when  considering  the  expressions  of  ET  and  ES  when  the 
service  times  are  all  exponentially  distributed. 

Remark  II.3.1.  Assume  the  servers  to  all  have  exponential  service  distributions. 
In  that  case 

Tk{x)  =  i-e-^k~Xpk)x  and  Wk(x)  =  1  -  —  e-^k~Xpk)x  ,  l<k<K, 

Pk 

for  all  x  >  0,  and  equations  (2.3.4a,  c)  become 


ET  =  V  — — — r  (2.3.5a) 

Hk  ~  *Pk 


(2.3.5c) 


where 

lk  :=  :  \I\  =  k}  . 

Formula  (2.3.5c)  is  better  suited  for  computations  for  values  of  K  less  than  ap¬ 
proximately  15.  For  larger  values  of  K ,  numerical  integration  of  (2.3.56)  appears 
to  be  more  economical  and  accurate. 

Remark  II. 3. 2.  Comparison  of  the  formulas  (2.3.36)  and  (2.3.3c)  shows  that 
the  system  time  has  more  pleasing  properties  than  the  resequencing  delay.  This 
appears  to  be  a  general  phenomenon  in  the  study  of  resequencing  systems  (see 
[BGP]  and  [Var.a]).  For  instance,  in  contrast  to  (2.3.3c),  it  was  not  possible  to 
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eliminate  the  2V s  in  (2.3.36)  by  using  Lemma  II. 3. 2.  Similarly,  the  easiest  way  to 
obtain  ER  is  through  equation  (2.3.46). 

Remark  II. 3. 3.  The  probability  of  being  a  star  job,  i.e.,  of  not  being  rese¬ 

quenced,  is 

K  ,oo  K 

P(.)  =  ij(0)  =  zw  n  Wi(t)dTk(t), 

k-1  ■'°  1=1 

l*k 

and  has  no  simple  expression  in  general  although  special  cases  are  given  in  Chapter 
IV. 

Remark  II.3.4.  If 

W  :=  m&x{W1,...,WK}  , 

then  (2.3.3c)  and  (2.3.4c)  can  be  rewritten  as 

1  dW 

S(x)  =  W{x)  -  -—(a;)  i  *  >  0  , 
whenever  W  is  differentiable,  and 

ES  =  EW  +  y  P(system  not  empty)  . 

Thus,  the  equation  for  S(-)  is  similar  to  (2.3.2)  and  seems  to  suggest  that  the 
system  admits  an  equivalent  M/GI/1  representation,  in  the  sense  that  S  =  W  +  B 
for  some  RV  B  independent  of  W.  Unfortunately,  such  a  representation  does  not 
seem  possible. 

Remark  II. 3. 5.  The  expressions  for  T,  R  and  S  depend  only  on  the  waiting 
and  response  time  distributions  of  the  parallel  M/GI/1  queueing  systems,  which 
are  well  known  at  least  by  their  Laplace  transforms.  If  the  distributions  Rjt(-) 
are  of  phase  type  (PH-type),  then  so  are  the  distributions  TVjt(-)  and  T*(-)  (see 
Appendix  II).  Since  the  maximum  of  independent  PH-type  distributions  is  again 
of  PH-type  [Neu.b],  the  distributions  of  T,  R  and  S  will  also  be  of  PH-type,  in 
view  of  Theorem  II. 3.1.  Although,  the  dimensionality  problem  in  obtaining  the 
maximum  of  PH-type  RVs  limits  such  computations  to  small  values  of  K ,  Theorem 
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II. 3.1  may  serve  to  guide  approximations  in  obtaining  the  performance  measures 
of  interest  for  larger  values  of  I\ . 

II. 4.  OPTIMIZATION  OF  THE  SYSTEM 

The  problem  of  finding  the  optimal  vector  of  switching  probabilities  that  min¬ 
imizes  the  expected  response  and  system  times  is  studied  in  subsections  II. 4.1  and 
II. 4. 2,  respectively.  Comparison  of  both  optimization  problems  provides  a  better 
comprehension  of  the  effects  of  resequencing. 

II. 4.1.  Minimizing  the  Response  Time 

The  nonlinear  program 


min  ET(pu.-.,Pk) 

(pi,—  ,Pk)GV 

has  already  been  addressed  in  various  forms  in  the  literature,  e.g.,  the  Capacity 
Assignment  problem  in  [Kle.b]  or  [Kel,  Ch.4|.  In  the  present  context,  it  has  been 
solved  by  Buzen  and  Chen  [BuC].  Nevertheless,  it  is  briefly  discussed  here  in  order 
to  compare  its  solution  to  the  corresponding  minimization  problem  for  ES. 

The  following  theorem  states  that  the  RV  T  is  stochastically  convex  with 
respect  to  the  load  allocation  vector  in  the  strong  st  sense  defined  in  Appendix 
II,  thus  providing  a  rationale  to  the  intuitive  arguments  given  for  the  algorithm  in 
[BuC]. 

Theorem  II. 4.1.  For  all  x  >  0,  the  mapping  p  T(p,  x)  is  (strictly)  concave  on 
V,  he,  {T(p),  peT>}e  SCX(st).  ■ 

Proof.  It  suffices  to  show  that  for  every  x  >  0,  the  Hessian  matrix  of  the  mapping 
p  T(p,  a:)  =  PkTk(pk,x)  is  negative  definite  on  {p  G  [0, 1]K ;  A pk  <  pk,  1  < 

k  <  K}  D  V.  Note  that  this  Hessian  is  diagonal.  Furthermore,  the  diagonal 
elements  are  strictly  negative  since  by  Theorem  A.II.4  each  one  of  the  mappings 
p*.  i-»  Tk(pk,  x),  1  <  k  <  K ,  is  concave  and  decreasing  for  all  x  >  0,  and  the  result 
thus  follows. 
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□ 


The  following  corollary  is  now  immediate  by  Definition  A.II.l. 

Corollary  II. 4. 2.  If  A  <  p,  the  constrained  optimization  problem  min^gx)  ET(p) 
has  a  unique  solution. 

Denote  this  optimum  by  p*  =  (p*, . .  .,p*K),  and  let  y  be  the  unique  root  of 
the  nonincreasing  function 


K 


Pk 


/(*)  =  1"Ea 


fc=l 


1  +  a\p\ 


1  - 


'Ip.kX  -  1  +  a\p\ 


defined  for  x  >  maxi<fc<ft-[(l  —  °\p\)l 2/x*]+.  It  is  an  easy  exercise  to  show  that  if 
y  >  1/Pfc  for  all  1  <  k  <  K,  then  a  straightforward  Lagrange  analysis  shows  that 
p*  is  given  by 


Pk 


Pk 

A 


1  - 


i  + 


2pky  -  1  +  alul 


1  <  k  <  K, 


(2.4.1) 


and  that  it  lies  in  the  set  V. 

On  the  other  hand,  if  1/pk  <  y,  then  p£  given  in  (2.4.1)  is  negative,  and  by 
Theorem  II. 4.1  the  minimum  is  located  on  the  boundary  of  T>,  i.e.,  at  least  one 
of  the  pit’s  is  0.  The  dimension  of  the  problem  is  thus  reduced  by  at  least  one 
and  this  leads  to  the  algorithm  given  in  Appendix  I  to  compute  the  vector  p*. 
Several  monotonicity  results  that  yield  considerable  computational  savings  in  this 
algorithm  are  also  presented  in  Appendix  I. 

Remark  II. 4.1.  When  the  service  time  distributions  are  exponential,  <7k/jk  =  1 
and  (2.4.1)  reduces  to 

V.  =  W  -  Vw  J;k~  Xl-  ■  1  <h<K.  (2.4.2) 

This  dramatically  simplifies  the  computation  of  p*  as  indicated  in  Appendix  I. 
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Remark  II. 4. 2.  Another  simplification  occurs  when  all  the  servers  are  identical 
(but  otherwise  general),  in  which  case 


Pk  - 


JL_ 

K' 


1  <  k  <  K  , 


so  that  sharing  the  load  equally  among  the  servers  is  optimal. 

The  following  theorem  generalizes  this  last  remark. 

Theorem  II.4.3.  When  the  service  time  distributions  are  all  identical,  the  •prob¬ 
ability  vector  p*  =  (1/K, . . . ,  1/K)  stochastically  minimizes  the  RV  T(p),  i.e., 
T{p*)  <st  T(p)  for  every  p  in  V. 

Proof.  Since  the  service  distributions  are  all  identical,  the  distribution  functions 
Tjt(-)  differ  only  by  the  input  rate  \pk ,  i.e.,  Tk(pk-,x)  =  Ti(pk,x),  1  <  k  <  K. 
Therefore,  for  all  x  >0, 


K  K  K  1 

T(p,x)  =  X>r,(M,s)  <  r,(^p?,x)  <  T,(J2  jj.i)  =  T(p-,x)  . 


k= 1 


Jt=l 


k=l 


The  first  inequality  follows  from  the  concavity  of  the  mapping  p  *-»■  T\{p,x)  for 
every  x  >  0,  while  the  second  inequality  is  a  consequence  of  the  fact  that 


K 

argmin{^p|;  p  G  V) 
k= l 


since  the  IR+  i->  IR+  mapping  p  ->  Ti(p,x)  is  monotone  decreasing  by  Lemma 
A.II.4. 

□ 


II.4.2.  Minimizing  the  System  Time 

The  complexity  of  the  expression  (2.3.4c)  makes  the  solution  to  the  problem 

min  ES(pi,. . .  ,pk) 

{p\,-,pk)£T> 
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much  more  difficult.  Even  in  the  exponential  case,  the  problem  is  too  complicated 
to  be  solved  “by  hand”  except  for  the  case  K  =  2.  In  this  case,  a  solution  is 
presented  in  the  form  of  the  root  of  a  5th  order  polynomial.  However,  the  numerical 
results  obtained  from  the  solution  of  this  special  case  lead  to  the  idea  of  a  simple 
asymptotic  approximation  which  generalizes  to  the  case  K  >  2.  This  is  discussed 
in  Section  II. 5. 

Another  case  where  the  optimization  problem  can  be  solved  is  when  all  the 
servers  are  identical,  but  otherwise  general.  For  this  configuration,  it  is  shown  that 
equal  routing  probabilities  still  achieve  the  minimal  average  system  time. 

II.4.2a.  The  Case  of  Two  Exponential  Servers 

In  the  rest  of  this  chapter,  when  K  =  2,  pi  >  p2  is  assumed  without  loss 
of  generality,  and  the  notation  p  =  pi  and  p  =  1  —  p\  =  P2  is  used.  The  set  V 
can  then  be  parametrized  by  the  scalar  p  so  that  the  statement  up  €  X>”  is  now 
equivalent  to  “(p,p)  €  2?”. 

With  this  notation,  the  average  system  time  given  in  (2.3.5c)  takes  the  form 

1  +  1  App(Mi  +  M2) 

Mi  -A p  p2-  Ap  M2  Mi  MiM2(Mi  +  M2  -  A)  ' 

Differentiating  (2.4.3)  leads  after  some  simplifications  to  the  relation 

tpbs^  = 

with 

D(p)  =  MiM2(Mi  -  Ap)2(M2  -  ap)2(m  1  +  M2  -  A) 

and 

N(p)  =  (mi  -  A p)2(M2  -  A p)2  a(p)  +  A/ui^2(mi  +  M2  -  A)2  b(p)  , 

where  the  linear  increasing  functions  a(p)  and  b(jp)  are  given  respectively  by 
a{p)  =  (mi  +  M2  -  A)(M2  -  Mi)  “  A(1  -  2 p)(/n  +  //2) 

and 


N(P) 

D(p) 


(2.4.3) 
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b(p)  =  p2  ~  Pi  ~  A(1  -  2 p)  . 


Note  that  D(p )  >  0  under  the  stability  assumptions  (2.2.1). 
The  scalars  a0  and  bo  given  by 


a0 


bo 


Pi  ~  P‘2 
2(pi  +  P2) 


and 


Pi  -  p2 
2A 


are  the  unique  roots  of  the  equations  a(p)  =  0  and  b(p)  =  0,  respectively.  Since 
Pi  >  P2  and  A  <  p\  +  p2,  it  is  easy  to  check  that  |  <  ao  <  60  <  pi/X. 

Lemma  II. 4. 4.  The  fifth  order  numerator  polynomial  N(p)  has  a  unique  root 
r 0  in  the  interval  T  =  (1  —  ^2.,^-).  Furthermore  this  root  lies  in  the  interval 

[max{a0)l  -  } , fe0 ] - 

Proof.  First,  the  following  observations  are  made: 

(i)  For  all  p  <  ao,  a(p)  <  0  and  b(p)  <  0,  so  that  N(p)  <  0,  while  for  all  p  >  bo, 
N(p)  >  0  by  the  same  arguments,  with  N(p)  =  0  if  and  only  if  p  =  a0  —  b0. 

(ii)  For  p  <  b0 ,  (1  —  C  T,  since  b0  <  Pi/X- 

(iii)  The  derivative  N'(p )  is  given  by 

N'(p)  =  —2X(pi  -  \p)(p2  -  Xp)a(p)b(p)  +  a'(p)(pi  -  \pf(p2  -  A pf 

+  b'(p)\pip2(pi  +  P2  -  A)2  . 

Two  cases  have  to  be  studied: 

If  1  _  £2.  <  ao,  then  [ao,&o]  C  T  from  (ii).  According  to  (iii),  N'(p)  >  0  for 
every  p  in  (a0,  bo),  since  there  a(p)h(p)  <  0  and  the  last  two  terms  of  N'(p)  are 
always  positive.  The  result  then  follows  from  (i)  since  N(p)  is  monotone  increasing 
in  [a0,&o]- 

On  the  other  hand,  if  1  —  >  a0,  then  N'(p )  >  0  for  every  p  in  [1  —  60] 

by  a  similar  reasoning.  Therefore,  since  iV(l  —  ^2.)  <  0  and  N(p)  >  0  for  p  >  bo 
by  invoking  (i),  it  is  plain  that  N(p)  has  a  unique  root  r0  in  (1  —  60]  C  T. 

The  result  thus  follows  by  combining  these  two  cases. 
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Corollary  II.4.5.  If  p  >  the  optimization  problem 


min  ES(pup2) 

(PuP2)€V 


has  a  unique  solution  (p'npl),  where 


=  min-fljT'o}  and  p\  =  1  —  p\ 


Proof.  If  r0  >  1,  then  N(p)  <  0  for  all  p  in  V.  Therefore,  ES(p )  decreases 
monotonically  in  V,  and  is  minimal  when  p  —  1. 

On  the  other  hand,  if  rQ  <  1,  then  ES(p )  is  minimal  at  p  =  r0. 


Remark  II. 4. 3.  When  K  —  2  with  exponential  servers,  the  average  resequencing 
time  takes  the  form 

ER(p)  =  £  +  p  __£__£_  Xpp(p  1  + 

Pl  -  \p  p2~  A P  p2  Pl  PlP2(Pl  +  P2  -  A) 

Although  the  form  of  ER(p)  is  very  similar  to  that  of  ES(p )  given  in  (2.4.3),  the 
optimization  problem  for  ER(p)  is  quite  different.  For  instance,  when  p\  =  P2, 
elementary  arguments  show  that  p  =  1/2  is  a  local  minimum  (resp.  maximum)  for 
ER  when  p  >  y/2—  1  (resp.  p  <  \/2*- 1).  Furthermore,  for  p  >  1/2,  p  =  1/2  is  the 
only  solution  of  the  equation  d  ER(p)/dp  =  0,  whence  it  is  the  global  minimum. 
On  the  other  hand,  it  is  plain  that  for  p  <  1/2,  p  =  1  in  T>  is  the  global  minimum, 
since  0  =  E  12(1)  <  ER(  1/2).  Therefore, 


(1,  if  p  <  1/2, 
arg  min  E  R(p)  =  < 

1  1/2,  if  p  >  1/2  . 
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Remark  II. 4. 4.  When  p i  =  p 2,  ao  =  bo  =  1/2  =  p\  so  that  ET(p )  and  ES(p) 
are  both  minimized  when  the  load  is  equally  allocated  between  the  servers.  This 
was  not  evident  a  priori,  since  ER(p )  may  be  maximum  in  this  configuration,  e.g., 
as  would  be  the  case  when  p  <  \/2  —  1  (see  also  [JeM.b]  and  [IL.c]). 

II. 4. 2b.  Optimal  Routing  in  the  Homogeneous  Case 

In  this  subsection,  the  system  is  studied  when  all  the  servers  have  identically 
distributed  service  times.  It  has  been  noted  in  Remark  II. 4.1  that  making  all 
routing  probabilities  equal  minimize  the  average  response  time  ET  for  the  ho¬ 
mogeneous  system.  The  following  theorem  states  that  this  result  still  holds  for 
ES  and  generalizes  the  observation  made  in  Remark  II.4.4.  This  shows  that  the 
response  time  dominates  the  resequencing  time. 

Theorem  II.4.6.  When  the  service  time  distributions  are  all  identical,  the 
probability  vector  p*  =  (1/K, . . . ,  1/K)  minimizes  the  function  ES(p). 

Proof.  As  noted  in  the  proof  of  Theorem  II.4.3,  the  distribution  functions  W*(  ) 
differ  only  by  the  input  rate  \pk ,  so  that  Wk(pk ,  •)  =  Wi (Pk,-),  1  <  k  <  K. 
Equation  (2.3.4c)  can  be  rewritten  as 

ES(p)  =  /  (1  -  Wi(pk,x))dx  +  -(l  -  JJ(1  -  poPk))  ■  (2.4.4) 

*'°  k=l  k=  1 

It  is  shown  now  that  both  terms  in  the  right-hand  side  of  (2.4.4)  are  minimum  at 
p* ,  so  that  their  sum  is  therefore  minimal  at  this  point. 

Strict  convexity  of  the  function  x  log(l  —  pox)  and  Lemma  A.III.l  of 
Appendix  III  imply  that  the  function  log(nf=i(l  —  PoPk )),  whence  the  product 
nf=i(i  —  p0pk),  is  maximum  at  p*.  The  second  term  is  therefore  minimum  at  this 
point. 

On  the  other  hand,  it  is  easy  to  show  by  Lemma  A.II.4  of  Appendix  II  that  for 
every  x  >  0,  the  positive  function  p  1— ►  log  W\(p,  x)  is  strictly  concave.  Therefore, 
by  the  same  argument,  the  point  p*  maximizes  the  product  n£=i  W\(pk,x)  for  all 
x>0,  and  minimizes  the  integral  in  the  first  term. 
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Remark  II. 4. 5.  Since  p*  —  (-jb, . . . ,  maximizes  the  product  WT-(p*,  %) 

for  all  x  >  0,  the  RV  max{RT;  1  <  k  <  K }  is  stochastically  smallest  at  this  point, 
i.e.,  when  the  RVs  Wk ,  1  <  k  <  K,  are  i.i.d.  This  observation  and  (2.3.3c)  suggest 
that  a  stronger  (stochastic)  optimality  result  as  in  Theorem  II. 4. 3  also  holds  for 
the  RV  S(p).  This  is  currently  under  investigation. 

II. 4. 3.  Examples  and  Comparisons:  K  =  2 

This  section  aims  at  a  better  understanding  of  the  variation  of  the  optimum 
switching  probabilities  with  the  system  parameters.  For  the  case  K  =  2  with 
exponential  servers,  the  notation  p*(\,pi,p2)  and  pi(\,  pi, p2)  is  used  for  p*  and 
pi  given  in  (2.4.2)  and  Corollary  II.4.5,  respectively,  in  order  to  explicitly  indicate 
their  dependence  on  the  system  parameters.  Figure  II.  1  displays  the  sets 

r P  =  {(pi,p 2)  |  A p  <  pi,  A(1  -p)  <  p2  and  pJ(A,  pi,  p2)  =  p  } 


and 

CP  -  {(^1,^2)  |  A p<  pi,  A(1  —  p)  <  p2  and  p|(A, pi,p2)  =p  }  , 

for  0  <  p  <  1  in  the  (pi,p2)  plane  when  A  =  1. 

Only  the  sets  Tp  and  Cp,  for  p  ranging  from  0.5  to  1  with  increments  of  0.05, 
are  drawn  in  Figures  Il.la  and  II. lb,  respectively.  The  sets  for  values  of  p  smaller 
than  0.5  follow  from  symmetry  by  interchanging  pi  and  p2.  The  dashed  curves  are 
the  lines  p\  +  p2  =  A  (the  stability  limit),  pi  =  p2  (To. 5  and  C0.5),  p\  =  p2  —  2A 
(the  asymptote  for  the  boundary  curve  of  Ti)  and  p\  =  p2  —  A  (the  asymptote  for 
the  boundary  curve  of  C\). 

An  important  observation  in  Figures  Il.la  and  II. lb  is  that  the  sets  Tj,  and 
Cp,  0  <  p  <  1,  are  smooth  curves.  These  curves  start  from  the  point  (A p,  A(1  —p)), 
since  the  set  T>  shrinks  to  the  point  (pi,P2)  =  (pi/AbP2/p)  as  p  tends  to  A  (high 
utilization).  This  provides  a  trivial  approximation  for  a  heavily  loaded  system. 
The  points  located  under  the  lowest  curve  in  Figures  Il.la  and  II. lb  are  points 
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of  the  sets  Ti  and  C\,  respectively,  and  correspond  to  systems  where  the  slower 
server  is  not  used  at  all,  i.e.,  p\  —  1  or  p\  —  1. 

Figure  II.2  (resp.  II. 3)  provide  the  same  comparison  for  a  DS  where  an 
M/M/1  queue  is  in  parallel  with  an  M/D/1  (resp.  M/Geo/1)  queue.  The  abbre¬ 
viation  “Geo”  denotes  a  geometric  distribution,  defined  by  P(B  —  kA)  =  qk(l  —  q), 
where  A  is  a  time  constant  and  q  a  fixed  probability.  For  this  distribution 
fx  =  (1  —  q)/Aq  and  pPa2  =  (1  +  A//).  The  parameter  A  is  taken  to  be  1/A 
in  Figure  II. 3.  These  families  of  distributions  are  chosen  for  their  simple  relation¬ 
ships  between  a  and  /. i . 

The  sets  Yp  and  Cp,  for  p  ranging  from  0  to  1  with  increments  of  0.1  are  drawn 
in  Figures  II. 2  and  II. 3.  As  in  Figure  II.  1,  points  located  under  the  lowest  curve 
are  points  of  lb  (or  C\ )  and  correspond  to  systems  for  which  the  second  server  is 
not  used.  By  symmetry,  the  curve  on  the  top  is  the  boundary  of  the  region  To  (or 
Co),  and  all  points  located  above  it  correspond  to  systems  where  the  first  server 
is  not  used  at  the  optimum. 

The  points  of  Tp  are  obtained  by  solving  (2.4.16)  numerically.  The  curves  Cp 
are  obtained  as  follows:  Since  the  first  server  has  exponentially  distributed  service 
times,  Wi(-)  is  as  given  in  Remark  II. 3.1.  Therefore,  after  some  simplifications 
and  using  the  notation  of  Section  II.4.2a,  ES  reduces  to 


ES(p)  =  ET2(p)+^--^~ 

Pl  p  2 


App  +  A p  W2*(p i  -  A p) 

P1P2  Pi  Pi  -  A p 


where  W% (•)  is  the  Laplace  transform  of  W2(-).  The  curves  Cp  are  obtained  by 
solving  the  equation  d  ES(p)/dp  =  0  numerically. 

The  comparison  of  the  curves  and  Cp  in  Figures  II.1-II.3  reveals  that 
P*  £  Pi  for  all  values  of  the  system  parameters,  i.e.,  the  presence  of  resequencing 
always  decreases  the  optimal  utilization  of  the  slower  server. 
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Figure  Il.la 

The  curves  Tp  for  p  =  0.5, 0.55, . . . ,  0.95, 1.0 
M/M/l  in  parallel  with  M/M/ 1 
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Figure  Il.lb 

The  curves  Cf  for  p  =  0.5, 0.55, . . . ,  0.95, 1.0 
M/M/ 1  in  parallel  with  M/M/ 1 
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Figure  II.2a 

The  curves  T,  for  p  *=  0, 0.1, . . . ,  0.9, 1.0 
M/M/ 1  in  parallel  with  M/P/1 
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Figure  II.2b 

The  curves  Cf  for  p  *  0, 0.1, . . .  ,0.9, 1.0 
M/M/1  in  parallel  with  M/P/1 
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Figure  II.3a 

The  curves  for  p  =  0, 0.1, . . . ,  0.9, 1.0 

M/M/1  in  parallel  with  M/Geo/1 
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Figure  H.3b 

The  curve®  Cp  for  p  «  0, 0.1, ...» 0.9, 1.0 
M/M/ 1  in  parallel  with  M/Geo/ 1 
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II.5.  APPROXIMATIONS  FOR  EXPONENTIAL  SERVERS 


This  section  develops  approximations  for  the  optimal  allocation  probabilities 
pf  and  p*  when  all  the  service  times  are  exponentially  distributed.  The  approxi¬ 
mations  are  motivated  from  Figures  Il.la  and  ILlb,  where  the  level  curves  Fp  and 
Cp,  0  <  p  <  1,  have  asymptotes  that  are  parallel  to  the  line  p\  =  More  impor¬ 
tantly,  both  the  curves  Cp  and  their  asymptotes  remain  close  to  one  another  for  all 
values  of  the  system  parameters.  Therefore,  the  idea  behind  the  approximations 
is  to  replace  the  curves  Cp  by  their  asymptotes. 

The  approximations  for  the  case  K  —  2  are  obtained  as  a  natural  extension 
to  the  discussion  in  Section  II.4.2  and  are  illustrated  separately  as  they  provide  a 
better  understanding  of  the  key  ideas.  Extensions  to  the  case  of  K  >  2  parallel 
servers  are  then  provided  through  the  asymptotic  expansion  of  the  formula  (2.3.5c). 
Although  a  simple  algorithm  is  available  in  Appendix  I  for  computing  the  vector 
p* ,  a  similar  approximation  is  also  derived  in  order  to  see  the  effect  of  resequencing 
more  clearly.  The  accuracy  of  these  approximations  is  validated  for  a  wide  range 
of  system  parameters. 

II. 5.1.  The  Case  K  —  2 


In  this  subsection,  the  notation  and  definitions  given  in  Section  II. 4. 2a  are 
again  adopted.  The  asymptotes  to  the  curves  Cp  in  Figure  I. la  can  be  obtained 
directly  by  noting  that 

a(p)  =  Oi  +  P2 )Kp)  ~  Kp*  ~  Pi) 


and  rewritting  the  equation  defining  Cp  as 
N(p) 


P  i 


=  0  =  b(p) 


(i  +  «)(i- pA)* («  _  +  ti  _  A)*' 

Pi  Pi  Pi  Pi  Pi  Pi  Pi  Pi  . 

_A(£i_i)(i-pA)2(a_pA)2. 

Pi  Pi  Pi  Pi 

Letting  pi  go  to  +oo  with  pifpi  going  to  some  constant  a,  it  follows  that  a  =  1, 
so  that  the  line 

Kp)  =  P2-Pi-  Ml  -  2p)  =  0 
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is  the  asymptote  to  the  curve  Cp.  As  observed  from  Figure  II. la,  the  asymptote 
Kp)  —  0  is  parallel  to  the  pi  =  p2  line  and  passes  through  the  point  (A p,  A(1  —p)). 
Since  the  curve  Cp  also  passes  through  this  point,  Cp  approaches  b(p)  under  both 
heavy  and  light  loads,  explaining  why  the  approximations  are  good  for  all  values 
of  the  system  parameters  as  indicated  in  Table  II.  1. 

A  similar  argument  shows  that  the  line  p2— pi  — 2A(1— 2 p)  ~  0  is  the  asymptote 
to  the  curve  Tp.  The  approximations  p\  and  p\  to  p\  and  p\  are  therefore  given 

by 


and 


2  +  ^  Pi  ~  P2  <  A, 

1  otherwise, 


(2.5.1a) 


5+“^  if  li,- lit  <  2X, 
1  otherwise. 


(2.5.16) 


Remark  II. 5.1.  When  p\  =  p\,  both  queues  are  stable,  so  that  this  approxima¬ 
tion  may  be  used  for  all  values  of  pi  and  p2  satisfying  the  condition  A  <  p.  On 
the  other  hand,  the  approximation  p\  is  valid  only  if  2A  <  pi  +  3p2  • 

Remark  II. 5. 2.  Note  that  p\  =  bo,  whence  by  Lemma  II.4.4  this  is  always 

greater  than  p\.  On  the  other  hand,  p*  can  be  shown  to  be  smaller  than  p\. 

Remark  II. 5. 3.  The  approximations  are  exact  when  pi  =  p,2  and  get  better  as 
the  service  rates  get  closer  to  each  other. 

In  order  to  quantify  the  accuracy  of  the  approximations,  numerical  examples 
are  collected  in  Table  II.  1  by  setting  A  =  1.  ET(p^)  is  also  displayed  to  see  the 
effect  of  resequencing  on  the  performance  measure. 


-36- 


(hl,^) 

Pi 

Pi 

%e 

ES(Pf) 

ES(p<) 

%e 

ET(p<) 

(1,0.25) 

0.8693 

0.875 

0.66 

10.15 

10.19 

0.34 

7.77 

(1,0.50) 

0.7425 

0.750 

1.00 

5.12 

5.13 

0.15 

3.95 

(1,0.75) 

0.6196 

0.625 

0.87 

3.39 

3.40 

0.35 

2.66 

(1,1) 

0.5000 

0.500 

0.00 

2.50 

2.50 

0.00 

2.00 

(2.5, 1.5) 

0.9532 

1.000 

4.91 

0.665 

0.667 

0.31 

0.648 

(2.5,1.75) 

0.8405 

0.875 

4.10 

0.647 

0.648 

0.14 

0.607 

(2.5,2) 

0.7273 

0.750 

3.12 

0.619 

0.620 

0.05 

0.568 

(2.5,2.25) 

0.6138 

0.625 

1.83 

0.586 

0.586 

0.01 

0.532 

(2.5, 2.5) 

0.5000 

0.500 

0.00 

0.550 

0.550 

0.00 

0.500 

(10,9.00) 

0.9787 

1.000 

2.18 

0.1111 

0.1111 

6.0e-3 

0.1110 

(10,9.25) 

0.8592 

0.875 

1.84 

0.1108 

0.1108 

3.0e-3 

0.1095 

(10,9.50) 

0.7396 

0.750 

1.41 

0.1102 

0.1102 

1.0e-3 

0.1081 

(10,9.75) 

0.6198 

0.625 

0.83 

0.1092 

0.1092 

3.0e-4 

0.1067 

(10,10.0) 

0.5000 

0.500 

0.00 

0.1079 

0.1079 

0.0e-0 

0.1053 

Table  II.  1. 

Asymptotic  Approximations  for  K  =  2 

The  approximation  p\  yields  a  relative  error  of  less  than  5%.  Moreover,  the 
relative  errors  in  ES  are  much  smaller  than  1%.  For  p*,  this  approximation  is 
not  as  good  for  small  values  of  p\  and  pi,  for  which  the  curves  Tp  and  their 
asymptote  are  far  apart.  However-;  the  curves  in  Figures  II. 1-3  indicate  that 
such  approximations  for  p*  get  better  as  the  service  time  distributions  become 
less  variable.  Of  course,  the  easily  computable  closed  form  expression  given  by 
equation  (2.4.2)  can  always  be  used.  Comparison  of  ES(p^)  and  ET(p^)  shows 
that  resequencing  degrades  the  system  performance  considerably  for  p  >  0.5. 

II.5.2.  The  Case  K  >  2 
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For  A  fixed,  the  optimization  problems 


min^T(p)  and  min  ES(p) 
p£V  v  '  P ev  v  ' 


are  considered  in  this  section  when  each  pk  goes  to  infinity. 

It  is  an  easy  exercise  [Lue]  to  show  that  the  optimal  probability  vector 
P*{p i )  •  •  • ,  Pk)  (similarly  p\p  1 , . . . ,  pk))  satisfies  the  Lagrange  equations 

r\ 

~ET{p)~  7(1)  -7fc2)  =  0  ,  1  <  k  <  K 

and 

k= 1 

(2) 

where  the  constant  7^.  >  0,  1  <  k  <  K. 

Let  I*  :=  {k  :  p*k  >  0}.  Then,  7^  =  0  for  every  k  in  I*,  so  that  p*  satisfies 
the  equations 

kET{p) = kET(p)  (2-5-2) 

for  every  k  and  l  in  I*. 

Using  (2.3.5),  the  partial  derivatives  in  (2.5.2)  are  given  by 


d 

dpk 


ET(p) 


Pk 

(Pk  -  A pk)2 


1  <  k  <  K  . 


(2.5.3) 


For  p*(pi, . . .  ,Pk)  =  P>  equations  in  (2.5.2)  define  surfaces  in  the 
(pi , . . . ,  /^ic)-space.  As  in  the  previous  subsection,  the  idea  is  to  approximate  these 
surfaces  by  their  asymptotes  to  provide  asymptotically  exact  approximations  to 
p*.  For  the  purpose  of  obtaining  these  approximations,  let 

Pk  =  &kPl  +  8k  +  0(  ), 

Pi 

for  some  constants  ak  and  6k  with  0  <  a*  <  1,  1  <  k  <  K,  and  07  =  1  and 
61  =  0.  Also  assume,  without  any  loss  of  generality,  that  p\  >  p?  >  ...  >  pk , 
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so  that  p*  >  0,  since  pk  •>  Pi  implies  p l  >  p*  by  (2.4.2).  The  equations  for  the 
asymptotes  can  now  be  obtained  as  the  first  order  term  in  pi  in  the  asymptotic 
expansion  of  (2.5.2)  as  p\  — >  oo.  Using  (2.5.3)  in  (2.5.2),  it  is  readily  seen  that 

KPk-P*)  =  \p\ (oik-ai)  +  ^(Sk-Si)  +  o(l)  ,  kj'ml*  . 
Summing  this  equation  over  all  l  in  I*  leads  to 


A  Pl  = 


A 

\Im 


+  2h~ 


Eigj*  ^ 

2\I*\ 


1 

+  2m 


Oik  - 


E 


161- 


Oil 


1 1* 


+  o(l) 


k  in  I* 


after  further  rearrangements,  where  |/*|  is  the  cardinality  of  I*.  It  is  clear  from 
this  equation  that  unless  otk  =  E/e/*-0*/!^*!’  P*k  cannot  be  in  V  as  p\  goes  to 
infinity.  Therefore,  the  pi's  are  all  in  T>  if  and  only  if  all  the  at’s  for  k  in  I*  are 
equal  to  each  other.  Since  an  =  1,  ak  =  1  for  all  k  in  I*. 

When  Oik  =  1,  &k  —  Pk  —  P l  +  0(1/ p\)  and  an  easy  analysis  then  leads  to 


\  *  Pk  .  2 A  pi*  /  ■  r* 

Xpt  ~  T  +  2|/*|  +  o(1)  ’  “ 1  ’ 


(2.5.4) 


where  ///*  =  Ylkei*  Pk-  Therefore,  the  following  asymptotic  approximation  is 
proposed  for  p* 

lik  i  L  a  T* 

2X  +  K  1 


P* 


2X  ^  2  A|7* 

o  k  £  r  , 


(2.5.5) 


where  I*  {k  :  pi  >  0}. 

9- 

Since  the  set  I*  is  in  turn  determined  by  p* ,  the  approximation  (2.5.5)  defines 
an  implicit  equation  for  p*.  However,  p*  can  be  obtained  through  the  following 
simple  algorithm  with  c  =  2.  Algoritm  II.  1  makes  use  of  the  relations  p\  >  . . .  > 
p*K  by  (2.5.5)  under  the  assumed  order  on  pk,  1  <  k  <  K. 

Algorithm  II.  1. 

(i)  Set  l*-K 
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(ii)  Compute 


(iii)  If  (a)  pi  >  0:  Stop,  and  set 


Pk  = 


S  +  t(i-E!=iS)  k  =  i,...,t 


k  =  1  +  1,...,K 


(b)  pi  <  0:  Set  p*  =  0  and  /  —  and  go  to  (ii). 


Note  that  the  inequality 


Pk  > 


pi*  —  2A 


(2.5.6) 


necessarily  holds  for  all  k  in  I*.  As  in  the  case  K  =  2,  (2.5.6)  does  not  guarantee 
stability  of  the  system  and  the  approximation  (2.5.5)  for  p*  may  not  be  used  if  the 


additional  condition 


x<!k+lV%?>k) 


is  not  satisfied. 


For  the  resequencing  problem, 


■g-BSfp)  =  -i-  no  -  ^ 

dpk  Pk  Pj 


(2.5.7  a) 


•  €-T  \i& 


A pj  \  A  Api  +  YljeiiPj  ~  ^ Pj ) 


-L  Vr-iV+1  V  TT  2EL  yi  1 

k  k  M  «  «  GW  W  -  W)2 


=  --*E— +«(4)  • 

/ifc  frf  PkPi  Pi 


(2.5.76) 
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Therefore,  pi  satisfies  the  relations 


Hk 


AE 

i^k 


Pi 

VkPi 


+ 


o(4)  =  i-AV 


PIPi 


5 


which  rewrites  as 


A(pI  “  Pi)  =  Pk-Pl+^^  ~(P*  ~  A*0  +  °(!)  > 

& 

for  every  k  and  l  in  I*  :=  {k  :  >0}. 

By  an  argument  similar  to  the  one  given  above,  p],  are  all  in  V  only  if  =  1 
for  all  A;  in  I* ,  and  the  asymptotic  expansion  of  pi  is  then 

M  =  W  +  El'j'l  +  0(1)  ,  k  in  i»  .  (2.5.8) 

The  condition  for  p\  to  be  in  ii  now  reads  as 

fit  >  4fT p  ’  l<k<K  ■  (2.5.9) 

This  region  has  the  same  shape  as  the  one  defined  by  (2.5.6),  but  is  “twice  nar¬ 
rower”,  due  to  the  factor  2  in  (2.5.6).  This  is  a  direct  generalization  of  what  has 
been  obtained  for  the  case  K  =  2.  The  constraint  (2.5.9)  is  sufficient  for  the  vector 
p^  given  in  (2.5.8)  to  be  in  V  and  the  implicit  equations  defining  the  approximation 

r' 

pi  can  be  given  by 


£*.  _l  x.~?tL  k  g  Jt 
0  k  £  ft  , 


(2.5.10) 


where  I*  :=  {k  :  pi  >  0}.  The  approximation  pi  can  again  be  obtained  by 
Algorithm  II.  1  with  c  =  1. 
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Equations  (2.5.6)  and  (2.5.9)  indicate  that  p*k  and  p\  obtained  from  (2.5.4)  and 
(2.5.8),  respectively,  are  strictly  positive  only  if  the  corresponding  pk  is  “large  enough”. 
In  other  words,  p*k  (resp.  pk )  is  set  to  zero,  if  pk  does  not  satisfy  (2.5.6)  (resp.  (2.5.9)), 
and  (2.5.4)  (resp.  (2.5.8))  has  to  be  applied  to  the  reduced  set  of  queues.  Intuitively,  this 
procedure  reflects  the  fact  that,  when  the  load  is  sufficiently  low,  it  is  more  advantageous 
to  process  all  the  jobs  on  a  few  faster  servers.  Comparison  of  (2.5.6)  and  (2.5.10)  also 
shows  that  the  resequencing  requirement  tends  to  further  decrease  the  number  of  slower 
servers  used. 

The  randomly  selected  numerical  examples  in  Table  II.2  indicate  the  accuracy  of 
this  asymptotic  approximation  for  I(  =  3  and  A  =  1.  The  exact  solution  is  obtained 
by  solving  the  corresponding  nonlinear  programming  problem  numerically.  As  in  the 
case  K  —  2,  the  approximation  ES(p *)  seems  to  be  an  upper  bound  to  ES(p *),  with  a 
maximum  relative  error  of  only  about  1%. 


(Pl,P2,P3) 

V 1 

ES(p t) 

ES(p f) 

%e 

(1.0, 0.8, 0.5) 

(0.549,0.360,0.091) 

(0.567,0.367,0.067) 

3.027 

3.039 

0.40 

(1.0, 0.5, 0.1) 

(0.743,0.257,0.000) 

(0.750,0.250,0.000) 

5.118 

5.125 

0.14 

(1.0, 0.2, 0.2) 

(0.855,0.073,0.073) 

(0.867,0.067,0.067) 

10.01 

10.10 

0.90 

(1.0, 0.1, 0.1) 

(0.926,0.037,0.037) 

(0.933,0.033,0.033) 

19.77 

19.99 

1.11 

(0.7, 0.4, 0.3) 

(0.562,0.267,0.171) 

(0.567,0.267,0.167) 

11.85 

11.86 

0.08 

(0.7, 0.3, 0.2) 

(0.631,0.233,0.135) 

(0.633,0.233,0.133) 

24.86 

24.89 

0.12 

(0.7, .25, .15) 

(0.666,0.217,0.117) 

(0.667,0.217,0.117) 

51.70 

51.72 

0.04 

(1.2, 1.0, 0.3) 

(0.594,0.406,0.000) 

(0.600,0.400,0.000) 

2.032 

2.033 

0.05 

(2.0, 1.0, 1.0) 

(0.938,0.031,0.031) 

(4.000,0.000,0.000) 

0.992 

1.000 

0.81 

(2.0, 2.0, 1.0) 

(0.500,0.500,0.000) 

(0.500,0.500,0.000) 

0.750 

0.750 

_ 

0.00 

Table  II.2. 


Asymptotic  Approximations  for  K  =  3 


CHAPTER  III 


QUASI-STATIC  LOAD  ALLOCATION  IN 
PARALLEL  QUEUES  WITH  RESEQUENCING 

III.l.  INTRODUCTION 

The  model  of  Chapter  II  is  again  considered  when  the  service  times  are  ex¬ 
ponentially  distributed  with  rate  /Xfc  in  queue  k ,  1  <  k  <  I\.  The  Poisson  arrival 
rate  is  still  denoted  by  A  and  the  notation  =  /ifc/A,  1  <  k  <  K,  is  used. 

In  this  chapter,  the  system  parameters  bk,  1  <  k  <  A”,  are  all  assumed 
unknown.  The  approximate  optimal  probability  vector  p*  of  the  Bernoulli  switch 
that  minimizes  the  average  system  delay  is  given  in  Section  II. 5  as  the  projection 
Vp  (p)  of  the  vector  p  given  by 

p  —  b  +  -jrr(l  —  bT e)  e  (3.1.1) 

onto  the  probability  simplex  V,  where  b  =  (6j , . . . ,  &#).  The  simple  form  of  (3.1.1) 
is  utilized  here  in  a  stochastic  approximation  (SA)  algorithm  for  computing  the 
vector  ph 

The  next  section  provides  a  brief  introduction  to  SAs  and  a  framework  for  the 
proposed  algorithm.  In  Section  3,  an  idle  time  measurement  model  is  described  to 
estimate  the  system  parameters  (see  also  [Bok]).  The  proposed  SA  algorithm  is 
given  in  Section  4  and  several  numerical  examples  are  collected  in  Section  5.  All 
the  RVs  in  this  chapter  are  defined  on  a  probability  space  (O,^-,  V). 
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III. 2  BACKGROUND  ON  STOCHASTIC  APPROXIMATIONS 

Let  h  be  a  function  from  IRK  — >  IR^,  and  consider  the  problem  of  finding  the 
unique  vector  p*  such  that 

h(p*)  =  0  • 

There  are  various  iterative  methods  for  determining  p*,  such  as  the  Newton’s 
method  or  its  variants.  These  methods  require  the  evaluation  of  the  function  h(-) 
and  its  derivatives.  In  many  situations,  unfortunately  one  can  only  observe  h(p) 
corrupted  with  error  or  noise,  i.e.,  only  h{p)  +  £  is  available  for  some  RV  £. 

Robbins  and  Monro  [RoM]  suggested  the  following  SA  procedure  for  the  so¬ 
lution  of  this  problem:  Given  an  arbitrary  initial  point  po^  and  a  decreasing 
sequence  of  positive  numbers  {an  n  =  0, 1, . . .}  such  that 

oo  oo 

Y.  a2n  <  oo  and  Y,  an  =  oo  ,  (3.2.1) 

n=0  n=0 

set 

Pn- f-1  —  Pn  i  (3.2.2) 

where  zn  is  the  outcome  of  the  measurement  at  the  ntfl  iterate  given  by 

zn  —  /i(pn)  +  t |n  ■>  n  =  0, 1, . . .  ,  (3.2.3) 

where  {£„,  n  =  0,1,...}  is  a  sequence  of  conditionally  independent  zero-mean 
RVs.  Robbins  and  Monro  proved  that  the  sequence  {pn,  n  =  0, 1, . . .}  obtained 
from  this  procedure  converge  to  p*  under  the  conditions  given  in  (3.2.1).  The  first 
condition  in  (3.2.1)  guarantees  that.the  jumps  pn+i  —  Pn  are  damped  to  achieve 
convergence,  while  the  second  condition  ensures  that  the  magnitude  of  the  jumps 
does  not  decrease  too  rapidly  to  allow  recovery  from  a  poor  choice  of  po . 

Since  [RoM],  SA  algorithms  became  increasingly  popular  in  applications  due 
to  their  ease  of  implementation  and  to  the  availability  of  a  comprehensive  theory 

In  this  chapter,  the  subscript  n  represents  the  iteration  number,  so  that  the 
kth  component  of  a  vector  xn  in  IR^  is  denoted  by  1  <  k  <  K. 
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concerning  their  asymptotic  behavior.  The  reader  in  referred  to  [KuC],  [HaN]  and 
the  references  therein. 

In  some  cases,  such  as  the  one  considered  here,  the  iterates  pn  have  to  lie  in 
some  constraint  set  E.  A  modified  version  of  the  basic  Robbins-Monro  algorithm 
is  given  by  Kushner  and  Clark  [KuC]  when  E  is  closed  and  bounded:  The  (n-f  l)st 
iterate  is  obtained  by  projecting  the  vector  pn  —  anzn  onto  the  set  E,  i.e., 

Pn+1  =  VsiPn  -  anzn)  .  (3.2.4) 

Sufficient  conditions  for  almost  sure  (a.s.)  convergence  of  this  modified  procedure 
to  p *  are  also  established  in  [KuC]. 

III. 3.  THE  MEASUREMENT  MODEL 

In  this  section,  a  simple  procedure  for  estimating  6*,  1  <  k  <  K,  from  idle 
time  measurements  in  the  kth  M/M/1  queue  is  described.  These  estimates  provide 
the  aforementioned  sequence  { zn ,  n  =  0,1,...}.  This  measurement  technique 
is  used  by  Kumar  [Kum]  and  Bonomi  and  Kumar  [BoK]  for  similar  quasi-static 
optimization  problems. 

Consider  a  stable  M/M/1  queue  with  utilization  p.  Assume  that  this  queue  is 
sampled  at  the  points  of  a  Poisson  process  with  rate  u,  independent  of  the  arrival 
and  service  processes  in  the  queue.  If  NT  is  the  number  of  samples  in  the  interval 
[0,  t],  and  rjT  is  the  number  of  times  the  queue  is  sampled  idle  in  this  interval,  then 
the  a.s.  relations 

lim  =  lim  —  =  1  —  p  (3.3.1) 

r— too  J\T  r— too  UT 

hold  (see  Appendix  IV). 

With  the  RV  yT  defined  by 

yT  :=  (1  -  2 ,  (3.3.2) 

UT 

it  is  plain  from  (3.3.1)  that 

lim  yT  —  p~l  ,  a.s. 
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Therefore,  yT  can  be  written  in  the  form 


yr  =  p  1  +  wT  (3.3.3) 

for  some  RV  wr,  where 

wT  — >  0  a. s.  as  r  -4  oo  .  (3.3.4) 

For  the  parallel  system,  assume  that  each  server  is  sampled  independently 
by  a  Poisson  process  with  rate  v.  After  every  measurement  interval  the  central 
scheduler  (Bernoulli  switch)  receives  the  number  of  times  each  server  was  found 
idle  by  the  Poisson  sampling  process.  Let  rjk,n,  1  <  k  <  K  and  n  =  0,1,..., 
denote  this  number  for  the  kth  server  in  the  nth  measurement  interval  of  length 
rn.  The  choice  of  the  deterministic  sequence  {rn,  n  =  0, 1, . . .}  is  crucial  for  the 
convergence  of  the  proposed  SA  algorithm  (see  Remark  III.4.3).  The  queue  length 
at  the  end  of  the  nth  measurement  interval  is  used  as  the  initial  queue  length  for 
the  n  +  lst  interval.  The  switching  probability  vector  pn  is  held  constant  during  the 
nth  measurement  interval  and  is  updated  upon  receipt  of  the  new  measurements. 
The  arrival  rate  into  queue  k  during  the  nth  interval  is  therefore  A Pk,n-  Define 

Vk,n  :=  ,  1  <k<K,  n  =  0,1,....  (3.3.5) 

VTn 

When  vrn  is  not  chosen  to  be  an  integer,  the  denominator  is  non-zero  and  yk,n  is 
always  well  defined. 

In  view  of  (3.3.2)  and  (3.3.3),  the  following  model  for  the  measurement  vector 
yn  with  components  yk,m  l  <  k  <  K,  is  proposed: 

Vk,n  =  bk  +  Wk,n  ,  1  <k<K,  n  =  0,1,...,  (3.3.6) 

where  the  noise  sequence  (i«n,  n  =  0, 1, . . .}  has  the  property 

wn  0  a.s.  as  rn  — >  oo  .  (3.3.7) 

III.4.  STOCHASTIC  APPROXIMATION  ALGORITHM 
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Let  the  function  h  :  Ut^  — >  IR^  be  defined  by 

1  K 

Kp)  :=p-b-  —  (1  -  'Yjbu)  e  ,  (3.4.1) 

k= 1 

and  define  the  sequence  of  vectors  { zn ,  n  =  0, 1, . . .}  by 


1  K 

n  :=  Pn  Un  ~  ^(1  ^  ^  l/fc,n)  C  ,  Tl  —  0, 1,  .  .  .  . 


(3.4.2) 


fc=i 


It  is  plain  from  (3.3.7)  and  (3.4.1)  that 

zn  —  h(jpn )  ti  0, 1, . . . 

where 


(3.4.3) 


If  h(p*)  =  0,  then  p*  =  V-pip*),  i-e.,  pt  is  the  projection  of  the  root  of  the 
continuous  function  h  onto  'D.  However,  since  the  vector  b  is  not  known  and  only 
yn  can  be  measured,  only  zn  =  h(pn)  +  (n  is  available.  This  setup  naturally  calls 
for  the  SA  algorithm  outlined  in  Section  2.  Since  the  vector  pn  needs  to  be  a 
probability  vector  in  order  to  obtain  the  measurements  for  the  (n  +  l)s<  interval, 
a  projection  algorithm  of  the  type  given  in  (3.2.4)  is  needed. 

To  summarize,  the  following  SA  algorithm  is  proposed  where  the  positive 
sequence  {an,  n  =  0, 1, . . .}  satisfies  the  conditions  in  (3.2.1): 

Algorithm  III.l. 


(i)  Start  with  a  probability  vector  po>  and  set  n  =  0. 

(ii)  Obtain  the  idle  time  measurements  1  <  k  <  K,  during  an  interval  of 

length  Tn. 

(iii)  Compute 


_  Pk,n 
Vk,n  —  ^  _  T)k,n  > 
VTn 


1  <  k  <  K  . 
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(iv)  Compute 


zn  ~  Pn  —  yn  —  JfO-  -  ^2  Vk'n )  e  • 

fc=l 

(v)  Compute  pn+1  -  Vu(pn  ~  anzn). 

(vi)  Set  n  <—  n  +  1,  and  go  to  (ii). 


Remark  III. 4.1.  Since  the  set  V  is  unknown,  in  step  (v)  the  projection  is  done 
onto  the  probability  simplex  U  :=  {p  £  [0, 1]^  :  Yl^=i  Pk  —  !}• 


Remark  III. 4. 2.  Note  from  (3.4.2)  that 


where 


bn 


—  (1  0,n)Pn  "h  anVn  • 


Therefore,  in  step  (v)  of  the  Algorithm  III.l,  pn+ 1  is  obtained  by  replacing  b  with 
bn  in  Algorithm  II.  1. 

Remark  III.4.3.  A  convergence  proof  for  the  Algorithm  III.l  is  currently  under 
investigation.  Numerical  observations  in  Section  5  indicate  that  the  algorithm 
converges  within  a  neighborhood  of  even  when  rn  =  r,  n  =  0, 1, . . ..  The  choice 
of  t  is  also  important  for  the  accuracy  and  the  robustness  of  the  algorithm.  Clearly, 
for  large  values  of  r  the  measurements  and  the  control  algorithm  derived  from  it 
will  be  more  accurate,  but  the  control  will  be  updated  less  frequently.  The  reader 
is  referred  to  [BoK]  for  a  discussion  on  choosing  r. 


The  sequence  {an,  n  =  0, 1, . . .}  is  typically  taken  to  be 

an  =  — 77 -  ,  n  =  0,l,...  (3.4.4) 

n  +  I 

for  some  a  >  0.  In  this  case,  the  asymptotic  normality  of  the  normalized  error  term 
^/n(pn  —pi)  has  been  established  in  the  literature  under  various  assumptions.  The 
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reader  may  find  an  account  of  the  relevant  references  in  [NeH].  Generally  speaking, 
if  {an,  n  =  0, 1, . . .}  is  of  the  form  (3.4.4),  then  limn_*oo  \fn  ( pn  —  pi)  is  normally 
distributed  with  zero  mean  and  covariance  matrix  C.  The  matrix  C  in  general 
depends  on  the  parameter  a  and  the  gradient  Vh(p*),  and  a  is  chosen  so  that  C  is 
minimal.  The  interested  reader  is  referred  to  [NeH,  p.  169]  for  the  general  form  of 
the  C  matrix.  For  the  function  h  given  in  (3.4.1),  V/i(p*)  is  the  identity  matrix, 
and  the  matrix  C  is  given,  after  some  elementary  computations,  by 

where  C$  is  the  covariance  matrix  of  the  RV  £oo-  Therefore,  C  is  minimal  when 
a  =  1,  and  the  sequence,  {un,  n  =  0, 1, . . .}  is  thus  taken  as 

an  —  .  1  j  ft  =  0, 1, ...  , 

n  +  1 

in  the  following  examples. 

III. 5.  NUMERICAL  EXAMPLES 

In  this  section  the  results  of  a  few  experiments  obtained  from  a  simulation 
program  that  implements  the  SA  Algorithm  III.l  with  rn  =  r,  n  =  0,1,...,  is 
presented.  The  numerical  examples  are  picked  among  the  ones  given  in  Table  II. 2 
of  Chapter  II,  i.e.,  K  =  3  and  A  =  1.  The  service  rates  Hk ,  k  =  1,2,3,  and  the 
corresponding  p*  vector  (from  Table  II. 2)  are  given  in  Table  III.l  for  each  example. 


Example 

P 

(P1,P2,P3) 

P T 

III.l 

0.25 

(2.0, 1.0, 1.0) 

(1.000,0.000,0.000) 

III. 2 

0.63 

(1.0, 0.5, 0.1) 

(0.750,0.250,0.000) 

III. 3 

0.91 

(0.7, .25, .15) 

(0.667,0.217,0.117) 

Table  III.l 

Test  Examples  for  the  SA  Algorithm 
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In  all  three  examples  r  =  1000  *  max{pk  1  :  k  =  1,2,3}  while  v  is  such  that 
vt  =  1000  (see  [BoK]).  In  all  cases  p0  =  (0.33,0.33,0.34).  The  curves  marked 
1,  2  and  3  in  Figures  III.  1-III.3  present  the  evolution  of  the  optimal  switching 
probabilities  pijn,  p2,n  and  p3)„,  for  the  first  10,  25  and  35  iterates,  respectively. 
Although  the  convergence  of  the  proposed  algorithm  is  not  established,  the  iterates 
converged  to  within  5%  of  the 'limiting  values  in  a  few  iterations. 
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•  100  -f-  1 


Figure  III.l. 

p„,  n  as  0, 1, ...  f  10  for  Example  III.l. 
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•  m  -r 1 

Figure  III.2. 

pn,  n  —  0,1,..., 25  for  Example  111.2. 
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Figure  III.3. 

pn>  n  =  0, 1, . . . ,  35  for  Example  III.3. 
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CHAPTER  IV 


ASYMPTOTIC  RESULTS  FOR  PARALLEL 
QUEUES  WITH  RESEQUENCING 

IV.l.  INTRODUCTION 

The  model  of  Chapter  II  is  again  considered  when  the  load  is  equally  allo¬ 
cated  to  K  parallel  queues  with  identical  service  time  distributions.  The  common 
service  time  distribution  is  denoted  by  B(-)  with  finite  mean  1/ji  so  that  the  sys¬ 
tem  capacity  is  now  K/j,.  When  the  load  is  equally  allocated  to  the  queues,  the 
distributions  of  the  waiting  and  response  times  are  identical  for  all  queues,  and  are 
denoted  by  W^(-)  and  respectively.  Note  in  this  case  that  the  RV 

is  also  the  disordering  delay.  The  system  time  S  and  the  resequencing  delay  R  are 
also  represented  by  and  R^K\  respectively,  to  indicate  their  dependence  on 
K. 

Attention  is  given  here  to  the  variation  of  the  RVs  T^K\  R^K>>  and  with 
K.  In  Section  2,  it  is  shown  that  the  RV  T'K>  is  stochastically  integer  convex 
and  decreasing  in  K,  while  ES^K^  is  decreasing  in  K  when  the  arrival  rate  to 
the  system  remains  constant.  Asymptotic  expressions  in  K  are  also  provided  for 
£<*>(.)  and  £(*>(•)  in  Section  2  to  establish  (i)  convergence  of  these 
distribution  functions  to  the  corresponding  distributions  in  the  M / GI /  oo  system 
with  resequencing  and  (ii)  asymptotic  monotonicity  and  integer  convexity  results 
for  these  RVs.  It  is  shown  that,  while  the  behavior  of  R^K^  in  general  depends  on 
the  load  of  the  system,  and  always  have  similar  structural  characteris¬ 
tics.  For  instance,  ES ^  is  also  (asymptotically)  integer  convex  and  decreasing 
in  K. 

In  Section  3,  the  arrival  rate  to  the  system  is  increased  linearly  with’  K.  A 
totally  different  limiting  behavior  is  now  observed  in  this  case:  R^  dominates 
T^k)  as  both  ER^  and  ES grow  as  log.R,  while  ET remains  constant. 
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IV.2.  ASYMPTOTIC  RESULTS  FOR  CONSTANT  LOAD 


In  this  section,  the  limiting  behavior  of  the  totally  homogeneous  system  is 
studied  when  the  arrival  rate  A  to  the  system  is  held  constant ,  and  the  number  K 
of  queues  is  increased.  The  RVs  and  correspond  to  the  waiting  and 

response  time  distributions  in  a  M/GI/1  system  with  arrival  rate  X/K  and  service 
time  distribution  B(-),  respectively.  The  notation  po  =  X/p  is  used  throughout 
this  chapter. 

The  following  result  follows  directly  from  the  convexity  result  in  Appendix  II 
and  extends  the  monotonicity  result  of  [Sto,  p.  82]. 

Corollary  IV.2.1.  In  the  homogeneous  system  with  equal  load  allocation, 
the  response  time  is  stochastically  integer  convex  and  decreasing  in  K,  i.e. 

{TW,  I<  €  IN}  €  SDCX(st). 

Remark  IV.2.1.  For  the  system  time  S ^ ,  only  asymptotic  results  are  currently 
available  (see  Corollaries  IV.2. 5  and  IV.2. 6). 

For  this  totally  homogeneous  system,  equations  (2.3.36)  and  (2.3.3c)  in  Chap¬ 
ter  II  can  be  rewritten  as 

pOO 

R(K)(x)  =  /  [W^K\x  +t)]K~1dTK(t)  (4.2.1  a) 

Jo 

and 

S^r)  =  [W^ix)]*  -  t4~[w{K\x)]K  (4.2.16) 

A  ax 


for  all  a:  >  0. 

As  K  goes  to  infinity,  the  load  in  each  queue  goes  to  zero,  and  the  distributions 
of  the  RVs  W('K')  and  tend  to  the  unit  step  function  and  to  B(-),  respectively. 
The  expansions  of  these  distributions  as  p  — *  0  is  the  subject  of  the  following 
lemma.  For  the  purpose  of  the  analysis  to  come,  the  expansion  of  the  waiting  time 
distribution  is  given  up  to  the  second  order  term,  while  a  first  order  expansion 
suffices  for  the  response  time  distribution. 
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Lemma  IV.2.2.  Let  Wp{-)  and  Tp(-)  be  the  distributions  of  the  waiting  and 
response  times  in  a  M/GI/1  queue  with  utilization  p  =  X/p,  where  X  is  the  arrival 
rate  and  1/p  is  the  mean  of  the  service  time  distribution  B(-).  Then ,  as  p  0, 

Wp(x)  =  1  -  p(l-V1(x))  +  p2(V2(x)-V \(x))  +  o(p2)  (4.2.2 a) 

and 

Tp(x )  =  B{x)  -  p(B(x)  -  V3(x))  +  o(p)  (4.2.2 b) 

for  all  x  >  0,  where 


(1  -  B(t))dt, 


V\(x  —  t)(  1  —  B(t))dt 


and 


V3(x)  =  p  f  B(x  —  <)(1  —  B(t))dt  . 

Jo 


(4.2.2c) 


Proof.  Note  that  Vi(-)  is  the  distribution  function  of  the  residual  service  time  Vi, 
and  that,  V2{-)  =  V\  *Vj(*)  with  *  denoting  the  convolution.  Equation  (4.2.2a)  thus 
easily  follows  by  inverting  the  first  two  terms  in  the  expansion  of  the  Pollaczek- 
Khinchin  transform  formula  [Kle.a,  p.  200] 

OO 

K(s)  =  (i-P)'£ptK(s)]k . 

k—Q 


and  by  rearranging  the  terms. 

Equation  (4.2.2b)  is  plain  since  Tp(-)  =  B  *  Wp{-)  and  V^-)  =  V\  *  B(-). 

D 

Remark  IV.2.2.  Note  that  in  a  M/GI/1  queue  in  statistical  equilibrium  a 
job  will  arrive  to  an  empty  system  with  probability  1  —  p,  and  that  there  will  be 
exactly  one  job  in  the  system  with  probability  p  +  o(p).  This  can  be  proved  using 
the  Pollaczek-Khinchin  formula  for  the  queue  length  distribution.  Therefore,  the 
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asymptotic  expansion  of  the  response  time  Tp  can  also  be  obtained  by  observing 


that 


TP 


'  B  with  probability  1  —  p, 

< 

V3  with  probability  p  +  o(p )  . 


The  following  asymptotic  expansions  for  the  distribution  functions  and  means 
of  the  RVs  T^k\  and  constitute  the  main  result  of  this  section. 

Theorem  IV. 2. 3.  Consider  the  system  of  K  homogeneous  -/GI/l  queues  in 
parallel  with  Poisson  arrivals  and  balanced  Bernoulli  loading  under  the  resequenc¬ 
ing  constraint.  As  I<  goes  to  infinity,  the  asymptotic  expansions  of  the  response, 
resequencing  and  system  time  distributions  are  given  by 


T^k\x)  =  B(x)  +  ^(V3(x)  -  B(x))  +  o(-i)  (4.2.3a) 

/•oo  _  1  fOO  — 

R(K\x)=  /  e^oVl(i:+t)dB(t)  +  ^  /  e-^x+t\F(x  +  t)-po)dB(t) 

Jo  K  Jo 

r  00  _  1 

+  ^  /  e-^x+t)dV3(t)  +  o(— )  (4.2.36) 


and 

S{K\x)  =  B{x)e~PoV^x)  +  ±e-poVl(x)  [(F(®)  -  Po)B(x)  +  p0V3(x)}  +  o{~ -) 

(4.2.3c) 


where 

2 

F(x)  =  pl(V2(x)  -  Vi(*))  +  poVi(x)  -  ^-Vi\x)  (4.2.3 d) 

for  all  x  >  0.  Their  means  are  given  by 


/•OO  _ 

ER(k)  =  /  (1  -  e~poV^x))B{x)dx  + 

Jo 


P  0 
2  K 


1  Po  —  (1  +  <?2P2) 


+  G 


(4.2.4a) 


+  C(-) 
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and 


(4.2.46) 


ES{k) 


1  -  e~p° 


+ 


i; 


(i 


,-po  Vi(x 


))dx  + 


2  K 


o  —  pO 


+  G 


+  °(^)  (4.2.4c) 


where 

G  =  J  (yl(x)  +  2(Vi(x)  -  V2(*)))  e~poV^x)dx  . 

,4 

Proof.  Since  denotes  the  response  time  in  a  M/GI / 1  system  with  arrival 
rate  X/K  and  service  time  distribution  B(-),  replacing  p  with  po/K  in  (4.2.26) 
gives  (4.2.3a).  Equation  (4.2.4a)  then  follows  by  simple  integration. 

Equations  (4.2.36,  c)  are  derived  from  (4.2.1a,  6)  and  (4.2.3a),  replacing 
W^K\x)  by  its  expansion  (4.2.2a),  with  A  replaced  by  X/K,  and  noting  that  as  K 
goes  to  infinity 


u  v  .  1  . 


K 


2v  —  u2 
2 1< 


+  4>)  . 


(4.2.5) 


Rather  than  using  (4.2.3c)  to  obtain  ES by  integration,  it  is  much  simpler 
to  start  from  the  relation 

ESrn  =  r  ( 1  -  [W^(*r)  dx  +  i  (l  -  (l  -  %f) 

derived  from  (2.4.4),  and  to  replace  W^K\x)  by  its  expansion  (4.2.3a).  Equation 
(4.2.4c)  then  follows  from  (4.2.5)  by  routine  manipulations.  Finally,  ER/k^  is 
again  computed  as  ES ^  —  ET^K\ 

□ 


Remark  IV. 2. 3.  The  terms  indicated  with  the  shorthand  notation “o(l/RT)”  in 
(4.2.3)  are  functions  of  x.  By  tedious  yet  straightforward  calculations  it  can  be 
shown  that  the  integrals  involving  these  functions  in  Theorem  IV. 2. 3  exist  and  are 
still  of  the  order  o(l/K),  so  that  the  expansions  given  in  the  theorem  are  valid. 
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Corollary  IV.2.4.  For  the  system  of  K  homogeneous  •/ GI/1  queues  in  parallel 
with  Poisson  arrivals  and  balanced  Bernoulli  loading  under  the  resequencing  con¬ 
straint,  the  limiting  distributions  of  the  response,  resequencing  and  system  times 
as  K  goes  to  infinity  exist  and  coincide  with  the  corresponding  distributions  of  the 
M/GI/o o  system  with  resequencing,  namely 


T°°(x)  =  B(x) 

(4.2.6a) 

R°°(x)  =  r  e~poV^x+tUB(t) 

Jo 

(4.2.66) 

and 

S°°(x)  =  B(x)e~poV^x) 

(4.2.6c) 

for  all  x 

>  0.  Their 

means  are  given  by 

ET°° 

1 

V 

(4.2.7a) 

ER°° 

/*00 

=  /  (1  -  e-poV^x))B{x)dx  =  ES°°  -  ET°° 

Jo 

(4.2.76) 

and 

ES°° 

=  1  “  *  **  +  /°°( 1  e-poVlix))dx  . 

*  Jo 

(4.2.7c) 

Proof.  Equations  (4.2.6)  and  (4.2.7)  easily  follow  by  letting  K  — *  oo  in  (4.2.3) 
and  (4.2.4),  respectively.  The  fact  that  these  distributions  coincide  with  the  cor¬ 
responding  quantities  in  the  M/GI/oo  system  with  resequencing  can  be  verified 
from  the  results  of  [HaP]. 

□ 

Remark  IV.2.4.  In  the  M/M/1  case,  the  asymptotic  expansions  (4.2.2a  —  b ) 
are  immediate  since 

Tp(x)  =  i-e~uG-P)x  =  1  -  e-^1  -  ptxxe-^  +  o(p) 

and 
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Wp (*)  =  1  -  pe-^-p)x  =  1  -  pe -  p2pxe~llx  +  o(p) 

for  all  x  >  0.  Therefore,  (4.2.4a  —  c)  simplify  as 

ETW  =  ET°°  +  +  o(4)  ,  (4.2.8a) 

K  p  K 

1  fPO  1  _  e-t  I 

ER{k)  =  ER°°  +  —  (1  -  e""0  -  2p0  +  2p0  ^  - - - dt)  +  o(—) 

(4.2.86) 


and 


rp0  1  -  e_< 


ESW  =  ES°°  +  — ^(1  -  e-'0  +  2p0  £  -f-dt)  +  o(j^)  .  (4.2.8c) 


When  the  service  times  are  deterministic,  (4.2.4a  —  c)  take  the  form 

etW  4 +  ^ +  4>  ■ 

£B(KI  =  A +  4> 


and 


^  \ + 


(4.2.9  a) 
(4.2.96) 

(4.2.9c) 


By  studying  the  sign  of  the  coefficient  functions  of  \/K  in  the  asymptotic  ex¬ 
pansions  of  Theorem  IV.2.3,  it  is  now  possible  to  give  asymptotic  integer  convexity 
and  monotonicity  results.  Additional  comments  are  also  made  in  Remark  IV.2.5 
without  the  tedious  details  in  calculations. 

Corollary  IV.2.5.  In  the  resequencing  system  of  K  homogeneous  -/GI/l  queues 
in  parallel  with  Poisson  arrivals  and  balanced  Bernoulli  loading,  the  mapping  K  ■— + 
ESW  is  asymptotically  integer  convex  and  decreasing. 

Proof.  The  mean  ES (/c)  is  asymptotically  integer  convex  and  decreasing  if  the 
term  in  square  brackets  in  (4.2.4c)  is  strictly  positive.  By  the  definitions  of  the 
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functions  Vi ( * ) ,  i  =  1,2,  given  in  (4.2.2c),  it  is  plain  that 

v2(x)=  /\(,-,W(t)  <  V?(x)  <  V^x) 

Jo 

for  all  x  >  0,  and  the  constant  G  is  thus  positive.  The  result  now  follows  easily 
since  e~p0/A  +  G  >  0. 

□ 

.j 

Remark  IV. 2. 5.  The  situation  for  is  more  complex  as  should  be  apparent 
from  (4.2.3b).  When  the  service  times  are  deterministic ,  it  is  plain  from  (4.2.9b) 
that  ERfK^  decrease  to  0  at  least  asymptotically.  In  the  exponential  case  however, 
the  study  of  the  term  in  parenthesis  in  (4.2.8b)  shows  that  asymptotically  ER .W 
increases  (resp.  decreases)  to  ER°°  for  po  <  pQ  (resp.  po  >  Po)  with  p£  ~  0.783652. 

Although  the  integer  convexity  and  monotonicity  of  ES^K^  is  given  only 
asymptotically  in  Corollary  IV. 2. 5,  the  monotonicity  of  ES ^  in  K  follows  from 
Theorem  II.4.6  for  all  K. 

Corollary  IV. 2.6.  In  the  totally  homogeneous  system,  the  expected  system  time 
ES decreases  with  K. 

Proof.  The  system  with  K  —  1  queues  can  be  obtained  from  the  system  with  K 
queues  by  setting  pk  =  1  /(K  —  1)  for  1  <  k  <  K  and  px  =  0.  The  result  therefore 
follows  from  Theorem  II. 4. 6  since 

esk  <  . jfrr.o)  =  esk- i , 

with  an  obvious  meaning  to  the  notation. 

□ 

Asymptotic  expansions  in  (4.2.3)  also  provide  the  following  weak  asymptotic 
stochastic  convexity  and  monotonicity  result  for  the  RV  S^K\ 

Corollary  IV. 2. 7.  If  is  the  end-to-end  delay  in  the  resequencing  system 
of  K  homogeneous  -/GI/l  queues  in  parallel  with  Poisson  arrivals  and  balanced 
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Bernoulli  loading,  then  for  all  x  >  0,  there  exists  a  finite  integer  K(x)  such  that  the 
mapping  (K(x),  oo)  i— ►  [0,1]  :  K  i— ►  >  x]  is  integer  convex  and  decreasing. 

Note  that  the  constant  K(x)  depends  on  x  since  the  term  o(\/K)  in  (4.2.3c)  also 
depend  on  i  as  indicated  in  Remark  IV.2.3. 

Proof.  The  result  follows  if  the  term  in  the  square  brackets  in  (4.2.3c)  is  strictly 
negative  for  all  x  >  0,  i.e.,  if  for  all  x  >  0,  K  S^K\x)  is  increasing  and  concave 
for  K  >  K(x).  By  using  the  definition  (4.2.3d)  of  F(- ),  this  term  can  be  rewritten 
as 

( F(x )  -  po)  B(x)  +  p0V3(x)  =p20  (V2(:r)  -  -  0  B(x) 

+  po(V3(x)-B(x)V1(x))  .  (4.2.10) 

It  was  shown  in  the  proof  of  Corollary  IV.2.5  that  V2  (x)  <  Vi(x).  Indeed,  it 
is  a  simple  exercise  to  show  that  V^ix)  <  V^{x),  unless  B(-)  is  the  unit  step 
function.  Therefore,  the  first  term  in  (4.2.10)  is  strictly  negative  for  all  x  >  0  since 
Vi(x)  <  1.  Similarly, 

y3(x)=  /  Bix-^dV^t)  <  Fi (x)B(x) 

Jo 

and  the  second  term  in  (4.2.10)  is  also  negative  for  all  x  >  0.  Therefore,  the 
coefficient  function  of  1/K  in  (4.2.3c)  is  strictly  negative  for  all  x  >  0,  and  the 
result  follows. 

□ 

Remark  IV. 2. 6.  The  stronger  asymptotic  result,  namely,  that  there  exits  a 
finite  integer  K*  such  that  the  mapping  (K*,o o)  ■-»  [0,1]  :  K  h*  P[S(k)  >  x]  is 
integer  convex  and  decreasing  for  all  x  >  0,  is  currently  under  investigation.  Note 
that  this  is  equivalent  to  the  RVs  {S^K\  K  >  K*}  being  SDCX(st). 

To  conclude  this  section,  asymptotic  expansions  of  the  probability  of  being  a 
star  job  are  given  in  some  special  cases.  In  the  totally  homogeneous  system,  the 
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formula  (4.2.1a)  takes  the  form 


Rw(0)  =  P(K)(*)  =  f  [W^ixyK-'clT^Xx). 

Jo 

In  the  exponential  case,  tedious  yet  elementary  computations  show  that 

Po  K 

and  the  asymptotic  result 


=—p  0 


Po 


2K 


pmv-lzfl  +  BgL+tL) 


thus  follows  from  (4.2.5).  Note  that  P°°(*)  =  (1  —  e~Po)/p0  can  also  be  derived 
from  the  results  of  [KKM]. 

In  the  deterministic  case,  the  asymptotic  result 


can  be  shown  to  hold.  Note  that  =  1  as  expected. 


IV. 3.  ASYMPTOTICS  FOR  INCREASING  LOAD 

In  this  section,  it  is  assumed  that  the  input  rate  Xk  into  the  homogeneous 
system  with  K  queues  is  A k  =  K A  for  some  fixed  A.  Thus,  the  load  of  each  queue 
remains  fixed  at  po  =  A/p  as  K  varies.  Therefore,  the  common  distributions  of 
the  i.i.d.  RVs  Wk  and  T*,  1  <  k  <  K,  are  independent  of  K  and  are  now  denoted 
by  W{-)  and  T(-),  respectively.  Equation  (2.3.4c)  now  reads  as 

ES{k>>  =  E(  max  Wk)  +  1  -  %~P^—  .  (4.3.1) 

xi<k<K  '  KX 


The  asymptotic  method  used  in  [BMS]  for  the  Fork-Join  queue  applies,  and  leads 
to  the  following  theorem. 
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Theorem  IV.3.1.  If  the  Laplace  transform  of  B(-)  is  rational ,  then 
ES{k)  =  7  log  I<(  1  +  o(  1))  =  ER(K)  +  ET  ,  K  oo 
-where  the  constant  7  depends  on  the  distribution  B(-). 

Proof.  Note  that  the  second  term  in  (4.3.1)  decreases  to  zero  with  rate  l/K. 
Under  the  enforced  assumption,  the  distribution  function  W(-)  also  has  a  rational 
Laplace  transform  and  the  well-knowil  asymptotic  result 

W(x)  =  1  -  Ce~qx{  1  +  o(l))  ,  a;  — >•  00  (4.3.2) 

thus  holds  for  some  positive  constants  C  and  q ;  see  [Bor,  p.  129]  for  a  similar 
proof  for  the  response  time  distribution.  Theorem  7.4  of  [BMS]  therefore  applies 
to  yield 

E(  max  W„)  =  ^(1  +  o(l))  , 

l<k<K  q 

and  the  theorem  follows  for  ES with  7  =  1  /q. 

The  second  equality  trivially  holds.  Note  that  since  ET  is  constant,  ER 
also  grow  as  log  K. 

□ 

Remark  IV.3.1.  This  limiting  behavior  of  the  resequencing  system  is  similar 
to  the  one  observed  in  a  Fork- Join  system  [BMS].  In  the  corresponding  Fork- Join 
system,  a  job  arriving  into  the  system  is  forked  into  K  tasks  and  each  task  is 
processed  in  one  of  the  K  parallel  queues.  As  soon  as  all  the  tasks  of  a  job  have 
been  serviced,  the  job  is  immediately  assembled  (i.e.,  tasks  are  joined)  and  leaves 
the  system.  When  the  service  times  are  i.i.d.  with  a  rational  Laplace  transform 
and  when  the  arrival  process  is  a  renewal  process,  the  moments  of  the  system  time 
in  this  Fork- Join  system  have  been  shown  to  grow  logarithmically  in  K  [BMS]. 

Remark  IV. 3. 2.  When  the  service  times  are  of  PH-type  with  representation 
(a,  A),  the  waiting  time  distribution  is  again  PH-type  with  representation  (poir,  L ), 


-  64  - 


where  i r  =  -/iccA-1  and  L  =  A  +  p0A°TT  [Neu.b].  In  that  case,  -q  in  (4.3.2)  is  the 
eigenvalue  with  largest  real  part  of  the  matrix  L  [Neu.b,  p.  63]. 

In  the  exponential  case,  L  =  \q  —  p,  so  that  7  =  l/(p  —  Ao),  and 

ESW  =  i^L(l  +  o(l))  .  (4.3.3) 

p  —  Xo 


Remark  IV. 3. 3.  When  the  load  in  each  queue  is  held  fixed,  the  limiting 

behavior  of  the  system  with  parallel  buffers  as  K  goes  to  infinity  differs  from  that 
of  the  M/GI/oo  queue  with  resequencing,  i.e.,  the  limiting  behavior  of  the  case 
when  there  is  a  common  buffer.  To  illustrate  this,  let  ES be  the  average  system 
time  when  the  DS  is  the  M/M/ 00  queue  with  arrival  rate  A K  and  service  rate  p. 
Replacing  po  by  Kpo  in  (4.2.7c)  then  yields 


ES 


OO 

K 


e  Kpoe  *  ^  dx  +  o(l) 

1  —  p~u 

- du  +  o(l) 

u 


I  (  f1 1 — Lldu  +  \og(Kp0)  -  / 
V  \J 0  «  J 1 


Kpo  e~u 


u 


du 


+  o(l) 


log  I< 


(l  +  o(l)) 


(4.3.4) 


where  the  last  equality  follows  since  ff°  e  u/u  du  <  00  so  that 


lim 

K—+  00 


log  p0K 


=  0  . 


Comparison  of  (4.3.3)  and  (4.3.4)  shows  that  although  both  ES^K^  and  ES^  grow 
logarithmically  in  K,  the  average  end-to-end  delay  has  a  smaller  growth  when  there 
is  a  common  buffer. 


Remark  IV. 3. 4.  At  increasing  load,  in  the  exponential  case,  the  probability  of 
being  a  star  job  is  given  by 


pl,°M  =  -  (!  -  P»)K] 


1  ,  1  X 

Kpo  +o(/i  )  ' 
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In  conclusion,  when  the  arrival  rate  into  the  system  is  also  increased  so  as 
to  keep  a  fix  load  to  each  queue,  the  delay  due  to  resequencing  grows  in  log  K 
while  the  disordering  delay  does  not  change  with  K.  This  result  indicates  that  the 
resequencing  delay  will  have  a  major  impact  on  the  end-to-end  delay  for  highly 
parallel  systems  with  reasonable  load  in  each  queue.  This  is  illustrated  in  Table 
IV.  1  when  the  service  time  distributions  are  exponential.  Table  IV.  1  presents  the 
ratio  ES^  jET  for  various  values  of  K  and  po.  As  apparent  from  this  table,  the 
resequencing  delay  is  much  larger  than  the  disordering  delay,  especially  for  large 
values  of  K  and  po. 


po\K 

10 

30 

50 

70 

ET 

0.1 

1.40 

1.99 

2.38 

2.67 

1.11 

0.2 

1.96 

2.87 

3.34 

3.66 

1.43 

0.5 

2.34 

3.34 

3.83 

4.15 

2.00 

0.7 

2.62 

3.65 

4.15 

4.48 

3.33 

0.9 

2.83 

3.89 

4.40 

5.24 

10.0 

Table  IV.  1. 


The  ratio  ESk/ET  for  exponential  servers 


-  66  - 


CHAPTER  V 


DYNAMIC  LOAD  ALLOCATION  IN  PARALLEL 
QUEUES  WITH  SYNCHRONIZATION 

V.l.  INTRODUCTION 

In  packet  switching  networks,  long  messages  are  broken  into  several  shorter 

.j 

packets  which  are  simultaneously  transmitted  from  source  to  destination  over  vir¬ 
tual  circuits  (parallel  channels).  Since  the  transmission  time  of  a  message  over 
a  channel  is  proportional  to  its  length,  this  “pipelining”  effect  can  considerably 
reduce  the  transmission  time  of  a  message  over  that  of  transmitting  the  message 
as  a  single  packet.  At  the  destination  node  the  messages  are  first  reassembled  and 
then  resequenced. 

This  chapter  considers  such  a  communication  system  when  the  parallel  chan¬ 
nels  are  identical.  The  DS  is  again  modeled  as  a  system  of  K  parallel  single  server 
queues  (channels).  The  assembly  and  resequencing  operations  of  jobs  (messages) 
are  both  performed  in  the  RB.  The  problem  of  dynamically  allocating  the  work¬ 
load  (message  length)  of  each  job  to  the  parallel  queues  is  considered  when  the 
interarrival  times  and  workloads  of  jobs  form  two  mutually  independent  sequences 
of  i.i.d.  RVs.  Specifically,  for  each  job  arriving  to  the  system,  the  problem  of  op¬ 
timally  partitioning  its  workload  into  L  <  K  tasks  (packets)  for  transmission  over 
the  parallel  queues  is  studied  when  the  cost-per-stage  is  taken  to  be  the  system 
time  of  a  job. 

The  problem  is  formally  defined  in  Section  2.  The  dynamic  programming 
methodology  is  used  in  Sections  3  and  4  to  obtain  the  optimal  allocation  policy 
which  minimizes  the  average  discounted  cost.  It  is  shown  that  the  optimal  policy 
is  Markov  stationary  and  steers  the  workload  in  each  queue  to  a  balanced  position 
as  fast  as  feasible.  In  Section  5,  the  optimal  policy  for  the  discounted  cost  is 
shown  to  also  minimize  the  corresponding  average  finite  horizon  and  the  long-run 
average  costs.  Section  6  considers  the  case  K  =  2  and  illustrates  how  the  solution 
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methodology  of  Sections  3-5  can  be  used  to  obtain  the  optimal  scheduling  policy 
when  the  messages  are  transmitted  as  a  single  packet.  Optimality  of  joining  the 
queue  with  the  smaller  workload  is  established  for  the  discounted,  finite  horizon 
and  the  long-run  average  costs. 

Denote  the  largest  component  of  any  vector  x  in  R^  by  ||a?||,  i.e.,  ||x||  := 
max{a:fc;  1  <  k  <  K).  A  function  /  :  R^  — *  R  is  said  to  be  monotone  in¬ 
creasing  if  x  <  y  implies  f(x)  <  f(y),  where  the  inequality  x  <  y  is  understood 

Af 

componentwise  as  Xk  <  yk,  1  <  k  <  K. 

V.2.  THE  MODEL  AND  THE  PROBLEM  FORMULATION 

The  model  considered  here  is  defined  on  some  probability  space  (SI,  IF,  P) 
that  carries  all  the  RVs  of  interest.  Jobs  arrive  to  a  system  of  K  parallel  identical 
queues  each  with  an  infinite  capacity  buffer.  The  time  between  the  nth  and  (n—l)3t 
arrivals  is  denoted  by  r(n),  n  =  1, 2, . . .  .  The  workload  of  the  nth  job  to  arrive  to 
the  system  is  represented  by  the  R+  -valued  RV  <r(n),  n  =  0, 1, . . .  ,  with  cr(0)  =  E. 
It  is  assumed  that  the  0th  job  arrives  to  the  system  at  time  t  —  0  and  finds  the 
initial  workload  in  the  system  to  be  the  R+ -valued  RV  W,  i.e.,  Wk  is  the  workload 
of  queue  k  at  t  =  0,  not  including  the  workload  of  the  0(/i  job. 

Upon  arrival,  a  job  is  partitioned  into  smaller  tasks  which  are  allocated  to  the 
parallel  queues.  Let  the  vector  u(n)  in  U  :=  {u  £  [0,1]^  :  uTe  =  1}  represent 
the  allocation  of  the  nth  job  to  the  queues,  i.e.,  cr(n)ujfc(n)  is  the  workload  of  the 
kth  task  of  the  nth  job,  1  <  k  <  K.  It  is  assumed  that  the  service  rate  (channel 
capacity)  in  each  queue  is  fixed  and  thus  equal  to  1  without  loss  of  generality. 
Therefore,  a(n)uk(n)  is  also  the  service  time  of  the  kth  task  of  the  nth  job.  After 
service  completion,  the  tasks  move  to  the  RB  to  await  service  completion  of  the 
other  tasks  that  belong  to  the  same  job.  After  all  the  tasks  of  a  job  complete  their 
service,  the  job  is  reassembled.  A  reassembled  job  further  awaits  in  the  RB  for  all 
the  jobs  that  have  arrived  to  the  system  earlier  to  be  reassembled. 

The  following  assumptions  are  made: 

(Al)  The  RVs  E  and  W,  and  the  sequences  of  RVs  {r(n),  n  —  1,2,. . .}  and 
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{<r(n),  n  =  0, 1, . . .}  are  all  mutually  independent. 

(A2)  The  RVs  {r(n),  n  =  1,2,...}  and  {cr(n),  n  —  0,1,...}  each  form  a 
sequence  of  i.i.d.  RVs  with  common  distribution  functions  A(-)  and  B(-), 
respectively,  with  the  properties  EA  <  oo  and  EB  <  oo. 

An  admissible  control  policy  i r  is  any  collection  7r  =  {7rn,  n  =  0,1,...}  of 
mappings 

7 rn  :  1R+  x  1R+  x(Wx  1R+  xE+)n->W.  n  =  0, 1, . . . 

The  collection  of  all  such  admissible  control  policies  is  denoted  in  the  sequel  by 
II.  For  every  7r  in  II,  the  1R+ -valued  RVs  {Wn(n),  n  =  0, 1, . . .}  and  the  W-valued 
RVs  {U*(n),  n  =  0, 1, . . .}  are  recursively  defined  by 

W*(n  +  1)  =  [W7r(n)  +  o{n)UK(n)  —  r(n  +  l)e  ]+  (5.2.1a) 

and 

U"(n  +  1)  =  7 rn+1  (E,  W;  U"(r ),  Ww(r  +  1),  *(r  +  1),  0  <  r  <  n)  ,  (5.2.16) 
with  initial  conditions 

W*(0)  =  W  and  17^(0)  =  7r0(W,E)  .  (5.2.1c) 

For  1  <  k  <  K  and  n  =  0, 1, . . .  ,  the  RV  W£ (n)  represents  the  workload  of  the 
kth  queue  at  the  nth  arrival  epoch,  while  the  RV  Uj?(n)  represents  the  fraction  of 
the  nth  job  allocated  to  the  kth  queue,  when  the  policy  7r  is  enforced. 

If  the  control  policy  tt  is  used,  then  the  sojourn  time  of  the  kth  task  of  the 
nth  job  is  Wj!{n)  +  U£(n)cr(n),  1  <  k  <  K,  and  the  system  time  57r(n)  of  the  nth 
job  is  thus  given  by 

S*{n)  =  || W*(n)  +  a(n)UK(n)\\  .  (5.2.2) 

Note  that  the  maximum  is  taken  over  all  queues,  even  if  some  of  the  components 
of  Un(n)  are  zero,  so  as  to  ensure  that  the  jobs  leave  the  system  in  sequence. 
Therefore  the  synchronization  delay  in  5w(n)  is  due  to  both  the  reassembly  and 
resequencing  operations. 
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For  every  /?,  0  <  0  <  1,  the  /3-discounted  cost  associated  with  a  policy  7 r  in 
II  is  defined  by 


:=  E 


,n= 0 


(5.2.3) 


and  the  minimization  problem  (P^)  of  interest  is  then 


(Pp) :  Minimize  Jp(ir)  over  II  . 

Let  A  and  B  be  two  independent  R+ -valued  RVs  defined  on  0,  with  probability 
distribution  functions  A(-)  and  P(-),  respectively.  By  simple  conditioning,  the 
relation 

P{W"(n  +  1)  <  w  ,  o(n  +  1)  <  o-  |  IF"}  (5.2.4) 

=  P{B  <  a}P  {[^"(n)  +  cr(n)I7"(n)  -  Ae]+  <  w} 

where 

IF"  =  a{W,  £}  V  cr{a(i),  Wn(i),  U”(i  -  1),  i  =  1, . . .  n} 

is  seen  to  hold  under  the  assumptions  (Al)  and  (A2),  for  every  admissible  policy 
7 r  in  II.  Therefore,  the  joint  distribution  of  the  pair  (W"(n  +  l),cr(n  +  1))  given 
IF"  depends  only  on  (W"(n),  <r(n))  and  U*(n).  This  suggests  that  the  problem 
(Pp)  can  be  viewed  as  a  Markov  decision  problem  with  augmented  state  process 
( W*(ri),cr(n )).  To  that  effect,  for  every  value  /?,  0  <  (3  <  1,  the  discounted 
cost-to-go  Jp  associated  with  an  arbitrary  policy  7T  in  II  is  defined  by 


Jg(w,cr)  =  E 


n=0 


W 


W, 


S  = 


(5.2.5) 


for  all  w  in  1R+  and  a  in  1R+.  The  corresponding  value  function  Vp  is  then  given 
by 


Vp(w,o )  =  mf  Jp(w,a)  . 
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V.3.  STRUCTURE  OF  THE  VALUE  FUNCTION 

For  any  mapping  /  :  1R+  x  JR+  i— »•  IR+,  introduce  Tftf,  u  in  U ,  and  Tpf  as 
the  mappings  1R+  x  BrL_)_  t— ►  1R+  defined  respectively  by 

Tpf(w,cr)  :=  ||id  +  au\\  +  f3E  [f  ([it;  +  au  -  Ae\+,B)]  (5.3.1) 

and 

Tpf(w,(r)  i±mmT$f(w,(r)  (5.3.2) 

for  all  w  in  IR+  and  tr  in  IR4..  The  expectation  in  (5.3.1)  is  taken  over  the  joint 
distribution  of  A  and  B. 

Theorem  V.3.1.  If  the  mapping  (w,<r)  i— ►  f(w,a)  is  monotone  increasing,  and 
the  mapping  w  t-+  f(w,a)  is  convex  for  every  a  in  1R+,  then  so  is  Tpf. 

Proof. 

Monotonicity:  Let  id1  and  w2  be  two  vectors  in  1R+,  and  a1  and  a2  be  two 
scalars  in  IR+.  If  wl  <  w 2  and  a1  <  cr2,  then  id1  +  <  it;2  +  &2u  for  every  u  in 

U,  and  the  inequalities 

Hit;1  +<71tt||  <  ||id2  +  <r2it||  and  [id1 +<71u  —  re]+  <  [w2  +  cr2  u  —  t  e]*  (5.3.3) 

readily  follow  for  every  r  in  IR+.  By  the  monotonicity  of  /,  the  second  inequality 
in  (5.3.3)  implies 

/([id1  +  ct1u  —  re]+,  cr')  <  /([id2  +  cr2u  —  re]+,  a')  (5.3.4) 

for  every  (r,  cr')  in  1R+  x  H+,  and  the  monotonicity  of  Tjff(w,  a )  now  follows  from 
the  first  inequality  in  (5.3.3),  i.e., 

t2)  (5.3.5) 

holds  for  every  u  in  U.  The  monotonicity  of  Tpf  is  now  immediate  from  (5.3.5). 
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Convexity:  Let  a  be  a  scalar  in  [0, 1].  First  the  mapping  (w,  u)  i— ►  f(w,  a)  will 
be  shown  to  be  convex  for  every  cr  in  1R+.  Fix  a  in  ILf,  and  let  u1  and  u2  (resp. 
w1  and  w2)  be  two  vectors  in  U  (resp.  IR.+ ).  The  relations 

|jcr(Q!U1  +  au2)  +  aw1  +  aw2\  |  =  Haleru1  +  to1)  +  a(cru2  +  io2)|| 

<  aWau1  +  to1 1|  +  a\\au2  +  io2||  (5.3.6) 

and 

[ctttf1  +  aw2  +  a(au 1  +  au2)  —  re]+  —  [a^u?1  +  cru1  —  re)  +  a(w2  +  cru 2  —  re)]+ 

<  afu;1  +  &U1  ~  re]+  +  a[w2  +  au2  —  re]+ 

(5.3.7) 

hold  for  every  r  in  It+,  owing  to  the  convexity  of  the  mappings  x  ||x||  and 
x  i— >  [x]+.  Therefore, 

/([aw1  +  aw2  +  ^(au1  +  au2)  —  re]+,  a') 

<  /(alw1  +  au1  -  re]+  +  a[w2  +  au2  -  re]+,  a') 

<  af([w 1  +  cru 1  —  re]+,cr')  +  af([w2  +  au2  —  re]+,a')  (5.3.8) 

hold  for  every  (r,  a1)  in  IR.-).  x  3EL_|_ .  The  first  inequality  in  (5.3.8)  follows  from  (5.3.7) 
and  the  monotonicity  of  /,  while  the  second  inequality  expresses  the  assumed 
convexity  of  /. 

The  convexity  of  Tfif  in  the  variable  (w,u)  thus  follows  from  (5.3.6)  and 
(5.3.8),  i.e.,  for  every  a  in  H+, 

T««l+««2  /(c^1  +  aw2 ,  a)  <  aTf  /( w 1 ,  a)  +  aTf  f(w2 ,  cr)  ,  (5.3.9) 

and  the  convexity  of  Tpf  in  w  now  follows  from  Lemma  A. III. 2  of  Appendix  III. 

□ 

The  key  optimality  result  for  problem  (Pp)  is  now  discussed. 

Theorem  V.3.2.  The  value  function  Vp  satisfies  the  Dynamic  Programming  equa¬ 
tion 

Vp  =  Tp  Vp  (5.3.10) 
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and  is  monotone  increasing ,  while  the  mapping  w  i— v  Vp(w,a)  is  convex  for  every 
a  in  JR+. 

Proof.  Let  {Vn,  n  =  —1,0,1,...}  be  a  sequence  of  mappings  x  11+  ]R+ 

defined  recursively  by 

Vn+i  =  TpVn,  n  =  —1,0,1,... 

where  V-i  is  defined  as  the  zero  mapping  on  H+  x  1R+.  Owing  to  the  non¬ 
negativity  of  the  cost  function  for  every  7r  in  II,  the  sequence  {Vn,  n  =  — 1,0,1,...} 
is  monotone  increasing,  and  the  convergence 

Voo (to,  cr)  :=  lim  Vn(w,cr)  (5.3.11) 

n—+  00 

thus  takes  place  for  all  w  in  IR+  and  a  in  IR+.  Note  from  the  continuity  of  the 
mapping  u  f-4  TpVn(w,cr)  that  the  level  sets  {u  6  U  :  Tp  Vn(w,a)  <  A}  are 
compact  for  every  w  in  ]R^,  a  in  3R+  and  A  in  IR.  Therefore,  it  follows  from 
the  recursive  definition  of  {Vn,  n  =  —1,0, 1, . . .}  and  the  Monotone  Convergence 
Theorem  that  Voo  —  Tp  V^.  Proposition  13  of  [Ber,  pp.  264-266]  thus  implies 
Vp  —  Voo  and  (5.3.10)  is  obtained.  The  second  part  of  the  proposition  is  now  an 
immediate  consequence  of  the  fact  that  monotonicity  and  convexity  are  preserved 
under  the  limiting  operation  in  (5.3.11). 

□ 

The  following  properties  of  the  value  function  is  an  immediate  consequence  of 
Theorem  V.3.2. 

Corollary  V.3.3.  The  value  function  Vp  exhibits  the  properties 

min  Vp(w,o')  =  V^(0, 0)  ,  (5.3.12a) 

(iu,<t)£R^  XR+ 

min  Vp(w,a)  —  Vp(0,cr)  ,  for  every  a  in  1R+,  (5.3.126) 

and 

min  Vp(w,o )  =  Vp(w,  0)  ,  for  every  w  in  1R+  .  (5.3.12c) 

<r€R  + 
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With  a  slight  abuse  of  notation,  let  tt*  be  any  Markov  stationary  policy  in  II, 
induced  by  the  mapping  7r*  :  1R+  x  R+  U,  which  satisfies  the  relation 

Vp(w,<r)  =  (rfM  Vp)  (w,a)  (5.3.13) 

for  all  w  in  R^  and  a  in  IR+.  Owing  to  the  compactness  of  the  control  set  U, 
there  always  exists  such  a  policy  7r*,  and  as  well  known  [Ber,  Prop. 13,  p.  264],  it 
is  optimal  for  problem  (Pp). 

The  following  property  of  the  value  function  will  be  crucial  in  identifying  the 
form  of  the  optimal  control. 

Lemma  V.3.4.  For  each  k,  1  <  k  <  K,  define  the  mappings  Rk  :  IR^  IR^  by 

{(wk+1,...,WK,wi,...,Wk)  1  <k<K 

(5.3.14) 

w  k=K  . 

With  this  notation ,  the  value  function  Vp  has  the  property 

Vp(w,  a)  =  Vp(Rk(w),  cr)  ,  1  <  k  <  K  ,  (5.3.15) 

for  every  ( w,a )  in  1R+  x  JR+ . 

In  fact  the  property  (5.3.15)  holds  for  any  permutation  of  w.  However,  the  cyclic 
rotations  Rk{w),  1  <  k  <  K,  suffice  for  our  purposes. 

Proof.  For  any  policy  7r  =  {7r„,  n  —  0, 1, . . .},  define  the  policies  Rk 7r,  l  <  k  <  K, 
by  RkK  =  {f?fc(7rn),  n  =  0,1,...}.  Let  na  be  the  set  of  Markov  stationary  policies 
in  n.  Since  the  queues  are  homogeneous  and  \\w\\  =  ||f?fc(i<;)||,  1  <  k  <  K,  it  is 
plain  that  the  relations 


JZ(w,a)  =  J*k*(Rk(w),a),  1  <  k  <  K,  (5.3.16) 

hold  for  every  n  in  ns,  w  in  1R^  and  a  in  1R+.  The  property  (5.3.15)  thus  follows 
by  noting  that  RkTl3  :=  {Rk tt,  7r  G  ns}  =  ns. 
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□ 


V.4.  THE  FORM  OF  THE  OPTIMAL  CONTROL 

In  order  to  provide  an  explicit  form  for  the  optimal  policy  7r*  given  by 


■K*(w,cr)  =  argmin  Tp  Vp(w,cr)  , 


for  all  w  in  R+  and  a  in  R+,  consider  first  the  vectors  u*(w,  a)  which  satisfy 


u*(w ,  a)  =  arg  min  Tp  Vp(tv,  a) 


(5.4.1) 


where  U'  :=  {u  €  R^  :  uTe  =  1}.  The  following  Lemma  can  be  proved  by 
arguments  very  similar  to  the  ones  given  in  Theorems  V.3.1  and  V.3.2. 

Theorem  V.4.1.  For  all  w  in  R+  and  a  in  R+,  the  TR,+  -valued  function 
u  TuVp(w,cr)  is  convex  on  the  convex  set  U' . 

The  optimal  policy  -ir*(w,a)  can  be  obtained  from  u*(w,a )  by  projecting  it 
onto  U,  owing  to  the  convexity  result  of  Theorem  V.4.1.  The  following  Lemma 
will  prove  useful  in  establishing  Theorem  V.4.3. 

Lemma  V.4.2.  If  the  mapping  tp  :  R^  R  is  convex  and  has  the  property  that 
tp(w)  =  tp(Rk(w))  for  all  1  <  k  <  K,  then 


min 

{wfzRK:  wT  e=c) 


tp{w)  =  tp(  — 


e) 


for  every  c  in  R. 

Proof.  For  every  vector  a  in  U  and  all  w  in  R^,  the  inequality 

ip  f  ^3  <XkRk(w)  \  <  J3  akHRk(w))  =  Hw)  (5.4.2) 

Vfe= 1  /  k= 1 

follows  from  the  assumptions  enforced  on  tp.  The  result  is  now  obtained  by  choosing 
a  =  e  in  (5.4.2)  since  j ?  i  Rk(w)  —  7^  (wre)  e  f°r  every  w  in  RA  . 


-  75 


□ 

Theorem  V.4.3.  For  every  a  in  11+  and  w  in  ]R^,  u*(w,cr )  in  (5.4-1)  is  given 

by 

u*(w,a)  =  ^  e~w  ■  (5.4.3) 

In  words,  u*(w,cr)  steers  the  workload  vector  to  a  balanced  configuration,  i.e., 
u*(w,  a)  is  such  that  cu*  +  w  has  e^ual  components.  Note  that  this  is  always 
possible  since  u*(w ,  a)  is  in  IR^. 

Proof.  Define  the  mapping  xp  :  1R+  JR+  by 

*(x):=\\x\ \\  +  pE[Vp([x-Ae]+,B)]  .  (5.4.4) 

J3y  Lemma  V.3.4,  the  equalities 

</■(«»(*))  =  ||B»(*)II  +  ^  [V„  ([«»(*)  -  A  e]+,B)] 

=  M+^E[V>(fl4([x-Xel+),B)] 

=  </>(x) 

hold  for  every  1  <  k  <  K.  Furthermore,  x  i— >  ip(x)  is  convex  owing  to  the  convexity 
of  the  functions  x  i->  ||x||,  [x]+  and  Vp(x,  •).  Therefore,  ip  satisfies  the  assumptions 
of  Lemma  V.4.2. 

For  every  w  in  IR^  and  o  in  1R+,  set  x(u)  =  w  +  cru.  Then,  equation  (5.4.1) 
can  be  rewritten  as 

u*(w,cr)  =  arg  min  ip(x(u))  . 

u£W 

Since,  for  every  u  in  U' ,  x(u)T e  =  <x  +  wTe,  i.e.,  x(u)T e  does  not  depend  on  u,  it 
follows  from  Lemma  V.4.2  that  u*(w,cr )  is  given  by 

w  +  au*(w,  a)  =  — (cr  +  wTe)  e  , 

and  (5.4.3)  thus  follows. 
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□ 

The  following  result  provides  a  necessary  and  sufficient  condition  for  u*(w,  a) 
given  in  (5.4.3)  to  be  in  U. 

Lemma  V.4.4.  If 

A  :=  ]T(|M|  ~wk)  , 

k= 1 

then  u*(w,a )  given  by  (5.4-3)  is  in  the  set  U  if  and  only  if  A  <  a. 

Proof.  Note  from  (5.4.3)  that  u*( w,  a)  lies  in  the  set  U  if  and  only  if 


0  < 


(a  +  wT  e ) 

K 


|«)| 


(5.4.5) 


since  u*(w,a)Te  =  1.  It  is  plain  from  the  definition  of  A  that 


-^(<7  +  wT  e)  -  |H|  =  -^(a  -  ]Tj(|HI  -  wk))  = 
and  the  result  holds. 

□ 

In  general,  when  u*(w,cr)  is  not  in  U,  the  convexity  result  in  Lemma  V.4.1. 
yields  the  following  result. 

Corollary  V.4.5.  For  every  a  in  H+  and  w  in  R+ ,  the  optimal  policy  for 
problem  (Pp)  is  given  by 


■K*(w,a)  =  Vu(u*(w,a)) 


(5.4.6) 


where  u*(w,cr)  is  given  in  (5.^.3). 

In  words,  the  optimum  allocation  strategy  7r*  steers  the  system  into  a  balanced 
configuration  as  quickly  as  feasible. 

The  following  corollary  is  now  immediate  and  states  that,  if  the  workloads  of 
all  the  queues  are  the  same,  then  the  optimum  allocation  strategy  tr*  keeps  the 
system  in  this  balanced  configuration. 
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Corollary  V.4.6.  For  every  a  in  IR-j.  and  w  in  such  that  w  =  c  e  for  some 
c  in  H+,  the  optimal  policy  is  given  by  7T*(w,a)  =  -b  e. 


Finally,  for  a  given  w  in  1FL+  and  a  in  R+,  an  algorithm  is  provided  for 
computing  the  vector  cr).  As  mentioned  before,  if  u*(w,cr)  is  not  in  U ,  then 

owing  to  the  convexity  of  the  mapping  U'  H+  :  u  TpVp(w,cr ),  7r*(w,a) 

is  on  the  boundary  of  U,  i.e.,  if  u*(w,a )  <  0  for  some  l,  then  7r*(w,a)  =  0. 
After  setting  a  component  of  tt*  to  zero,  the  problem  reduces  to  allocating  the 
workload  to  a  reduced  set  of  queues  so  that  the  remaining  components  of  tc* 
should  be  recalculated  from  (5.4.3)  and  (5.4.6).  The  following  result  states  that  if 
several  components  of  u*(w,a)  are  negative,  then  the  corresponding  components 
in  7r*(w,cr)  can  all  be  set  to  zero  at  once,  thus  facilitating  the  computations. 

Lemma  V.4.7.  For  every  non-empty  subset  E  of  {1,...,AT}  with  cardinality 
\E\  define  the  vector  u(E)  in  ]RjBl  by 

uk(E )  :=  -(c(£)  -  wk)  ,  keE  with  c(E)  =  — .  (5.4.7) 

a  \jbj\ 


Let  l  and  k  be  different  elements  of  E.  If  ui(E )  <  0,  then  Uk(E)  <  0  implies  that 
uk(E\{l})<  0. 

Proof.  The  result  follows  from  the  following  routine  calculations 


c(E  \  {/})  = 


+  ’Li&Ewi 

\E\-1 


wi  +  v  +  T,ieEwi 

\E\-1  \E\~1  \E\(\E\  - 1) 


a+J2ieEWi 

\E\ 


<  wk 


where  the  first  and  second  inequalities  follow  from  c(E)  <  u>i  and  c(E)  <  wk, 
respectively. 

□ 

Lemma  V.4.7  leads  to  the  following  algorithm 


Algorithm  V.l. 


-  78  - 


(i)  Set  E  «-  {1,2 . 

(ii)  Compute  the  vector  u(E)  from  (5.4.7). 

(iii.a)  If  Uk(E)  >  0  for  all  k  in  E ,  then  STOP.  The  optimum  allocation  vector  is 
given  by 


ir*k(w, 


Uk(E) 

0 


k  in  E 

k  in  {1  ,...,K}\E  . 


(iii.b)  Else,  for  every  k  in  E  such  that  Uk(E)  <  0,  set  E  <—  E  \  {A:}  and  go  to  step 
(ii). 


V.5.  THE  FINITE  HORIZON  AND  LONG-RUN  AVERAGE  COSTS 


In  this  section,  the  finite  horizon  and  the  long-run  average  cost  problems  are 
briefly  discussed.  It  is  shown  that  in  both  cases  the  optimal  policy  is  the  one  given 
for  the  discounted  cost  problem. 


V.5.1.  The  Finite  Horizon  Problem 

For  any  policy  ir  in  II  and  n  =  0, 1, . . .,  the  n-stage  total  and  average  expected 
costs  are  defined  respectively  by 


K  K<t) 


E  !»"(••) +  <K0t''(0ll 


L*=o 


W  =  w,  E  =  <7 


(5.5.1) 


and 


Jn(w,a) 


1 

n  + 1 


(5.5.2) 


for  all  w  in  1R+  and  a  in  IFLf.  The  corresponding  value  functions  are  given  respec¬ 
tively  by 

Vn{w,o)  =  mi  JZ{w, a) 

7rEiI 

and 


Vn(tc,cr)  =  inf  Jni,wi cr)  =  —y^Vn(w,v)  ■ 

Tren  n  +  i 

For  the  total  expected  cost  the  dynamic  programming  equation  is  given  by 
Vo (w,  a)  =  min  ||  w  +  au ||  (5.5.3) 

u£U 

Vm(w,  a)  =  min  {\\w  +  cru||  +  E  [Vm-i  ([ro  +  cru  -  Ae]+,  B)]  }  .  m  =  1, . . .  ,n 

u€l4 
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Let  7r^  be  the  vector  in  U  which  minimizes  the  right  hand  side  of  the  mth 
equation  in  (5.5.3)  for  m  =  0,1,..., n,  and  let  ir*  =  {7r^,  m  =  0,1,..., n)  be 
the  corresponding  optimum  policy.  The  following  lemma  shows  that  the  functions 
w  t— >  Vm(w,  a),  m  =  0, 1, . . . ,  n  all  enjoy  the  structural  properties  of  the  function 
Vp(w,  a).  Consequently,  an  argument  similar  to  the  one  given  in  Section  V.4  shows 
that  n*(w,  a)  given  in  (5.4.6)  is  optimum  for  the  (total  expected  cost)  finite  horizon 
problem  as  well.  To  that  end,  let  Rw  be  any  permutation  of  the  vector  w  in  1R^, 
and  let  i?”1  be  the  inverse  of  R,  i.e.,  R^^Rw)  =  R(R  1u>)  =  w. 

Lemma  V.5.1.  For  every  a  in  1R+,  ihe  mappings  w  t->  Vm(w,  cr),  m  =  0, 1, . . . ,  n 
are  convex  and  have  the  property  Vm(w,a)  =  Vm(Rw,a). 

The  following  Lemma  will  be  useful  in  proving  Lemma  V.5.1,  and  follows  from 
Lemma  A. III. 2  of  Appendix  III  by  observing  that  R(U)  =  U. 

Lemma  V.5.2.  If  the  mapping  rj>  :  x  U  1R+  is  jointly  convex  and 

has  the  property  ij>(Rw,Ru )  =  ij}(w,u)  for  all  permutation  R,  then  the  mapping 
<j)  :  • — ►  1R+  given  by 

<f>(w)  :=  min  ij>(w,u) 
u£U 

is  also  convex  and  has  the  property  =  <f>(Rw). 

Proof  of  Lemma  V.5.1.  The  result  follows  by  induction.  The  function  (w,  u ) 
ipo(w,  u)  =  \\w  +  (tu ||  clearly  satisfies  the  assumptions  of  Lemma  V.5.2.  The  result 
is  therefore  true  for  the  function  w  i— >  Vq(w,  a)  for  every  a  in  IR4-. 

For  m  =  1,2,...,  write 


Vm(iu,cr)  =  +  au ) 

u  £U 


where 

ipm{x)  =  ||a;||  +  E  [Vm_i  ([ x  -  Ae]+,B)]  . 

For  every  a  in  1R+,  if  the  function  Vm_1(-,cr)  is  convex  and  has  the  property 
Vm-i(R  ■ ,0 )  =  Vm_i(-  ,cr),  then  the  function  <j>m(w,u )  :=  +  au)  clearly 
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satisfies  the  conditions  of  Lemma  V.5.2.  The  function  w  Vm(w,  a)  therefore  has 
the  desired  properties,  thus  completing  the  induction  argument. 

□ 

Since  the  optimal  policy  tv*  for  the  total  expected  cost  is  independent  of  n,  it 
is  also  optimum  for  the  average  finite  horizon  problem.  This  can  also  be  seen  by 
writing  the  corresponding  dynamic  programming  equation  for  Vn  and  repeating 
the  arguments  presented  above. 

V.5.2.  The  Long-run  Average  Cost  Problem 

The  long-run  average  cost  for  any  policy  n  in  II  is  given  by 

J*(w,(t)  =  lim  J*(w,cr)  (5.5.4) 

n — MX) 

for  all  w  in  R^f  and  a  in  R4..  Since  J*(w,cr)  is  minimized-by  tv*  for  each  n, 
the  long-run  average  cost  is  therefore  also  minimized  by  tv*  .  In  this  section  the 
optimum  cost  Jn*(w,(r)  is  obtained  under  the  stability  condition  EB  <  K  EA. 
Let  the  subset  B  of  R^f  be  defined  by 

B  :=  {w  £  R^f  :  w  =  c  e  for  some  c  >  0}  . 

Also  define,  with  a  slight  abuse  of  notation,  the  work  process  {Wn(t),  t  >  0}  for 
every  tv  in  II,  where  W£(t)  is  the  workload  in  queue  k  at  time  t  under  the  policy 
l  <k  <K.  Note  that  with  this  definition,  if  tn  is  the  nth  arrival  epoch,  then 
W*(tn)  =  W*(n)  for  n  =  0, 1, . .  .. 

If  the  RV  T*  is  defined  by 

Tn  :=  inf{<  >  0  :  WK{t)  €  B)  , 


then 

Tn*  <  W  a.s.  (5.5.5) 

Note  that  TK  is  the  first  time  the  work  process  reaches  a  balanced  configuration 
under  the  policy  7v  in  II.  To  see  (5.5.5),  let  TJ  be  a  realization  of  T7r  for  to  — 
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(iy,cr(0),  r(l),cr(l), . . .)  in  fl.  If  ||u;|j  <  r(l),  then  then  the  system  will  be  empty 
at  time  ||u>||,  and  Tf  =  IMI-  On  the  other  hand,  since  7r*  steers  the  system  to  a 
balanced  position,  arrivals  into  the  system  will  only  result  the  system  to  reach  to 
a  balanced  position  earlier. 

Let  L-J  denote  the  integer  part  of  a  real  number.  The  sequence  of  EVs  {/n,  n  = 
0,1,...}  defined  by 

1  LT'*J 

fn  :==  ;  i  ^  S  (0  n  =  0, 1,... 

n  +  i  :=0 

is  uniformly  bounded  by  5Zl=o  by  virtue  of  (5.5.5)  and  converge  to  0  a.s. 

asn-4  oo.  It  is  plain  that  the  RV  Sn  (i)  has  finite  conditional  mean  given 

W.  Therefore,  for  some  c  and  o'  in  H+,  the  relation 

lim  J*  (ty,cr)=  lim  J *  (c  e,o')  (5.5.6) 

n— ^oo  n— ►oo 

holds  by  the  Bounded  Convergence  Theorem  [Bil,  p.  180]  for  every  w  in  1R^  and 
o  in  1R+. 

Under  the  assumed  stability  condition  EB  <  K  EA,  the  right  hand  side 
of  (5.5.6)  converges  to  the  average  response  time  in  the  GI/GI/1  queue  with 
interarrival  distribution  A( x)  and  the  service  time  distribution  B(Kx)  by  Corollary 
V.4.6.  The  following  theorem  summarizes  this  argument. 

Theorem  V.5.3.  The  policy  tt*  given  in  (5.4-6)  is  also  optimal  for  the  long-run 
average  problem.  Furthermore,  if  EB  <  K  EA,  then  the  corresponding  long-run 
average  cost  is  equal  to  the  average  response  time  in  the  GI/GI/1  queue  with 
interarrival  time  distribution  A(-)  and  the  service  time  distribution  B(K-). 

V.6.  OPTIMAL  SCHEDULING  WITH  NO  PIPELINING:  K  =  2 

This  section  considers  the  problem  formulated  in  Section  2  when  K  =  2  and 
when  the  control  set  W={«£  {0,  l}2  :  u\  +  u2  =  1},  be.,  when  the  incoming  jobs 
are  not  allowed  to  be  broken  into  smaller  tasks  and  are  scheduled  to  the  parallel 
servers  in  one  piece.  It  is  plain  that  the  structural  results  of  Section  3  still  hold  in 
this  case  and  that  the  optimal  policy  7r*  satisfying  (5.3.13)  is  Markov  stationary. 
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First  it  is  shown  that  the  vector  u*(w,cr )  in  W  =  {u  €  R2  :  «i  +  U2  =  1} 
given  by  (5.4.3)  can  also  be  used  to  characterize  the  optimal  scheduling  policy  in 
this  case.  For  every  e  in  R+,  let  u\  and  u2  be  the  two  vectors  in  U'  which  are  at 
a  distance  y/2  e  from  u*(w,cr),  i.e., 


and 


u)  =u*(w,  cr)  +  e 


(5.6.1a) 


=u*(w ,  cr)  +  e 


(5.6.16) 


The  following  result  states  that  the  function  U'  h->  1R-).  :  u  i-+  TuVp(w,cr)  is 
symmetric  around  u*(w,  a)  for  every  w  in  R^.  and  a  in  R+. 

Lemma  V.6.1.  If  the  vectors  u\  and  u2t  are  as  given  in  (5.6.1),  then 


T;‘V^w,a)=T^Vg(w,a) 


for  every  e  and  a  in  R+  and  w  in  R+. 

Proof.  Since  w  +  ou*(w ,  cr)  =  c  e  where  c  =  (a  +  wi  +  u>2)/2,  the  relations 

w  +  au\  =  c  e  +  ere  (  ^ 


and 


i  2  .  i  —  1 

w  +  crut  =  c  e  +  ere 


easily  follow  from  (5.6.1).  Therefore,  the  relation 


Ri(w  +  cru\ )  =  ili  (to  +  eru2) 


holds  for  the  permutation  ill  defined  in  (5.3.14).  The  result  now  follows  from 
Lemma  V.3.4.  and  the  definition  of  the  function  TpVp  in  (5.3.1). 

□ 
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It  is  clear  from  Lemma  V.6.1  that  it  is  optimal  to  schedule  the  arriving  job  to 
the  first  (resp.  second)  queue  if  u*(w ,  cr )  is  closer  to  ^  ^  ^  (resp.  ^  ^ )  than  ^  ^ 

(resp.  ^q^)-  The  following  result  is  therefore  obvious  given  the  form  of  u*(w,a ) 
in  (5.4.3). 

Theorem  V.6.2.  For  every  a  in  1R+  and  w  in  1R+,  the  optimal  policy  tt*  is  given 
by 


(0  < 


v*(w,a)  =  { 


10 


if  W\  <  lt> 2 


if  W\  >  W2  . 


In  words,  it  is  optimum  to  join  the  queue  with  the  smaller  workload. 

The  natural  extension  of  this  result  to  the  case  of  K  (>  2)  queues  is  currently 
under  investigation. 
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CONCLUSIONS  AND  FUTURE  RESEARCH 


In  this  thesis,  issues  of  performance  evaluation  and  optimal  routing  in  a  system 
with  K  parallel  queues  were  studied  under  the  resequencing  constraint.  In  Chap¬ 
ters  II-IV,  jobs  arrive  according  to  a  Poisson  process  and  are  allocated  randomly 
to  the  parallel  queues  by  a  Bernoulli  switch. 

In  Chapter  II,  the  service  time  distributions  at  different  queues  are  assumed  to 
be  independent  with  general  and  possibly  different  distributions.  The  system  under 
the  resequencing  constraint  is  compared  to  the  system  without  the  resequencing 
requirement.  It  is  established  that  the  presence  of  resequencing  severely  compli¬ 
cates  the  analysis  and  considerably  degrades  system  performance.  In  contrast  to 
the  simple  expressions  for  the  response  time  T,  the  expressions  for  the  distribu¬ 
tions  of  the  resequencing  and  system  delays  are  very  complicated.  In  particular, 
the  static  optimization  problem  of  choosing  the  Bernoulli  switching  probabilities 
to  minimize  the  average  system  time  ES  can  not  be  carried  out  by  hand  for  the 
general  case.  On  the  other  hand,  the  RV  T  is  shown  to  be  stochastically  con¬ 
vex  in  the  load  allocation  vector  in  the  st  sense  defined  in  Appendix  II,  and  the 
minimization  of  the  average  response  time  ET  is  carried  out  for  the  general  case. 

In  the  presence  of  resequencing,  two  special  cases  are  considered  for  the  op¬ 
timization  problem:  (i)  Identical  servers  and  (ii)  Exponential  servers.  When  the 
service  time  distributions  at  different  queues  are  identical,  ES  is  minimized  at 
equal  load  allocation,  and  at  this  optimal  configuration,  ES  is  decreasing  in  K. 
Stronger  results  are  obtained  for  the  RV  T.  Indeed,  equal  load  allocation  stochas¬ 
tically  minimize  T,  and  at  this  optimum  configuration,  T  is  integer  convex  and 
decreasing  in  K  in  the  st  sense.  The  weaker  result,  namely,  the  asymptotic  convex¬ 
ity  of  ES  in  I<  is  established  in  Chapter  IV.  Monotonicity  and  convexity  results 
in  the  st  sense  for  the  RV  S  are  currently  under  investigation  for  an  arbitrary  K. 
A  weaker  asymptotic  result  in  that  direction  is  included  in  Chapter  IV. 

For  the  exponential  service  time  distributions,  with  even  K  —  2,  the  optimum 
switching  probability  vector  does  not  admit  a  closed  form  expression  and  is  char- 
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acterized  by  the  unique  root  of  a  fifth  order  polynomial.  Nevertheless,  the  study 
of  the  optimization  problem  for  the  case  K  =  2  lead  to  simple  but  very  accurate 
approximations  for  K  >  2.  The  asymptotic  approximations  are  not  limited  to 
exponential  servers  and  can  be  applied  to  any  family  of  distributions  with  one 
parameter.  For  instance,  approximations  for  the  optimal  probabilities  displayed 
in  Figures  II.2  and  II. 3  may  be  obtained  by  the  same  technique.  Since  the  curves 
Cp  displayed  in  these  figures  are  closer  to  their  asymptotes  then  the  ones  drawn 
in  Figure  II.  1,  the  asymptotic  approximations  are  believed  to  be  even  better  for 
parallel  systems  with  less  variable  service  time  distributions. 

Both  the  approximate  and  the  exact  (when  applicable)  expressions  for  the  op¬ 
timum  switching  probabilities  indicate  that  the  faster  servers  have  to  carry  more 
traffic  both  with  and  without  the  resequencing  constraint,  but  that  the  resequenc¬ 
ing  requirement  tends  to  further  decrease  the  amount  of  the  traffic  switched  to  the 
slower  servers. 

In  Chapter  III,  the  simple  form  of  the  approximate  formula  is  used  in  a 
stochastic  approximation  algorithm  when  the  system  parameters  are  unknown. 
The  algorithm  makes  use  of  some  system  measurements  to  update  the  switching 
vector.  The  choice  of  the  lengths  of  the  measurement  intervals  crucially  affects  the 
rate  of  convergence  and  the  robustness  of  the  algorithm.  This  is  left  as  a  future 
research  topic. 

In  Chapter  IV,  the  asymptotic  behavior  of  various  system  delays  in  K  is 
studied  when  the  load  is  equally  allocated  to  the  homogeneous  servers.  When 
the  Poisson  arrival  rate  is  held  fixed,  asymptotic  expansions  are  provided  for  the 
distribution  functions  and  the  first  moments  of  the  RVs  T,  R  and  S.  These  asymp¬ 
totic  expressions  are  used  to  prove  convergence  of  the  system  statistics  to  those 
of  the  M/GI/oo  system  with  resequencing,  as  K  goes  to  infinity.  The  asymptotic 
expansions  are  also  used  to  establish  asymptotic  stochastic  monotonicity  and  con¬ 
vexity  of  the  RV  S  in  I<,  while  the  RV  R  asymptotically  may  increase  or  decrease 
depending  on  the  arrival  rate.  Therefore,  despite  the  different  behavior  of  R,  the 
RVs  T  and  S  have  (asymptotically)  similar  structural  characteristics. 
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In  the  case  when  the  arrival  rate  into  the  system  is  also  increased  so  as  to 
keep  a  fixed  load  to  each  queue,  the  resequencing  delay  dominates  the  queueing 
delay;  while  ET  remains  the  same,  ER  grows  to  infinity  in  log  K. 

In  Chapter  V,  the  parallel  system  is  considered  under  a  different  set  of  assump¬ 
tions  on  the  arrival  and  service  processes.  In  this  chapter,  jobs  arrive  according 
to  a  renewal  sequence  and  are  broken  into  smaller  tasks  for  processing  at  different 
queues.  The  servers  are  assumed  identical  with  fixed  capacities.  The  optimum, 
dynamic  load  allocation  problem  is  considered  when  the  workloads  of  the  jobs  are 
i.i.d.  and  independent  from  the  arrival  time  sequence.  The  cost  associated  for  allo¬ 
cating  each  job  is  chosen  to  be  its  end-to-end  delay,  including  the  synchronization 
delays  due  to  the  reassembly  and  resequencing  operations.  The  optimal  allocation 
policy  that  minimizes  both  the  long  run  average  and  the  discounted  costs  is  shown 
to  be  the  one  that  derives  the  workloads  in  the  queues  to  a  balanced  position  as 
fast  as  feasible.  The  same  optimization  problem  when  the  queues  have  different 
processing  capacities  is  currently  under  investigation. 

The  solution  to  the  optimization  problem  considered  in  Chapter  V  also  sheds 
light  into  the  characterization  of  the  optimum  scheduling  policy  when  the  jobs  are 
not  allowed  to  be  broken  into  smaller  tasks.  For  the  case  of  two  parallel  queues, 
scheduling  the  incoming  jobs  to  the  queue  with  the  smaller  workload  is  shown  to 
be  optimal.  The  optimality  of  this  policy  for  the  case  of  K  >  2  parallel  queues  is 
also  under  investigation. 
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APPENDIX  I 


COMPUTATION  OF  THE  OPTIMAL  PROBABILITY  VECTOR  p* 


In  this  Appendix  an  algorithm  for  computing  the  optimal  probability  vector 
p *  of  Section  II.4.1  is  given.  A  simplified  version  of  this  algorithm  is  also  provided 
when  the  service  time  distributions  are  exponential. 

As  noted  in  Section  II.4.1,  if  the  numbers  p£,  1  <  k  <  K,  obtained  from 
(2.4.1)  are  all  non-negative,  then  they  lie  in  the  set  V  and  constitute  the  solution 
to  the  optimization  problem  min {ET(p)  :  p  6  T>}.  On  the  other  hand,  if  one  or 
more  of  them  are  negative,  from  Theorem  II.4.1,  the  solution  is  on  the  boundary 
of  V ,  i.e.,  at  least  one  of  the  pk's  is  zero  and  equation  (2.4.1)  is  applied  to  the 
reduced  system.  The  following  argument  provides  a  computationally  efficient  way 
of  computing  p*  in  this  case. 

Let  the  functions  //,  1  <  l  <  K,  be  defined  by 

1  + 

2pkx  -  1  +  alul 

k*i 


for  x  >  maxi<jt<i<-[(l  -  er£p|)/2pfc]+.  The  functions  /  and  fh  1  <  l  <  K,  are  all 
decreasing  and  have  unique  zeros,  denoted  by  y  and  yi,  1  <  /  <  K,  respectively. 
Denote  p*k  in  equation  (2.4.1a)  by  p*k(y).  Note  that  p*k(y)  and  p*k(yi)  correspond 
to  the  optimal  probabilities  over  the  set  of  indices  E  =  {1, 2, . . . ,  K}  and  E  \  {/}, 
respectively.  The  following  lemma  leads  to  Algorithm  A.I.l  for  computing  the 
optimal  probability  vector  p*. 

Lemma  A.I.l.  Ifp*{y )  <  0  for  some  l  in  E,  then  pk(y)  <  0  for  some  k  in  E\{1} 
implies  pk{y\)  <  0. 

Proof.  Since  the  function  f\  is  decreasing  and  fi(y)  =  f(y )  +  P*(y)  =  P*(y)i  the 
condition  p*(y)  <  0  implies  yi  <  y.  Therefore  pk{yi )  <  Pkiv)  <  since  the 
functions  pk  are  all  increasing. 
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□ 

Algorithm  A.I.l. 

(i)  Set  E  <-  {1,2,..., if}. 

(ii)  Compute  y  as  the  solution  of  the  equation 

v"  Mfc  L  I  i  +  vl  1  =  1 

kts  X  [  V  2^kX  -  1  +  ak4  J 

(iii)  For  all  k  in  E,  compute  pi  as 


H  1  /  1  +  °llA 

\  j1  V  2m  -  1  +  aifi  _ 


Let  F  C  E  be  such  that  k  is  in  F  if  and  only  if  pi  <  0. 


(iv)  If  F  =  0,  then  STOP. 

Else,  set  E  *—  E\F,  pi  *—  0  for  every  k  in  F,  and  go  to  (ii). 


Note  that  the  monotonicity  of  the  functions  fi  can  be  used  to  computational  ad¬ 
vantage  in  step  (ii).  However,  this  step  is  still  the  most  computationally  intensive 
step  of  the  algorithm,  and  must  be  avoided  as  much  as  possible.  Lemma  A.I.l 
allows  for  the  index  set  E  to  be  reduced  to  E  \  F  in  step  (iv)  at  once,  instead  of 
reducing  E  one  element  at  a  time  by  stopping  when  a  negative  p*k  is  found.  This 
provides  an  improvement  over  the  algorithm  given  in  Buzen  and  Chen  [BuC].  Note 
that  pl(y)  >  0  does  not  imply  pk(yi)  >  0,  and  the  algorithm  needs  to  go  back  to 
step  (ii). 

In  the  exponential  case,  this  algorithm  can  be  improved  by  using  the  simplified 
formula  (2.4.2).  Assume,  without  loss  of  generality,  that  /^i  >  >  . . .  >  pk-  It 

is  easy  to  show  that  if  pi  >  0,  then  >  p*k ,  1  <  k  <  K.  Therefore,  if  one 

of  the  probabilities  is  positive,  so  are  the  ones  with  a  lower  index.  In  particular, 
if  P*k  ^  then  all  the  p*k  s  will  be  probabilities.  This  simplifies  step  (iii)  of 
Algorithm  A.I.l  and  leads  to  the  following  algorithm 
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Algorithm  A.I.2. 

(i)  Set  n  «—  K 
(ii)  Compute  p* : 


*  _  /fn  _  /  J2i= 1  Mi  _  -I  \  \/~Pn 

p”  A  1  A  Tr=i\^7 


(iii)  If  p*  <  0,  then  set  p*  <—  0,  n  «—  n  —  1,  and  go  to  (ii). 
If  P*n  >  0,  then  p*  =  (pj, . . .  ,p*,0, . . .  ,0),  where 


„*  Pk  /  ]Li=l  Pi  _  \ 

P‘  A  (  A  TLiVw’ 


1  <  jb  < 


n. 
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APPENDIX  II 


STOCHASTIC  CONVEXITY  RESULTS  FOR  THE  M/GI/1  QUEUE 

In  this  appendix,  a  notion  of  strong  stochastic  convexity  [SS.b]  is  defined  and  a 
stochastic  convexity  result  is  established  for  the  M/GI/1  queue  by  means  of  simple 
arguments.  This  convexity  result  states  that  the  stationary  waiting  and  response 
times  in  the  M/GI/1  queue  are  both  stochastically  increasing  (resp.  decreasing) 
and  convex  in  the  arrival  (resp.  service)  rate.  This  provides  an  important  step  in 
establishing  some  of  the  convexity  and  optimality  results  in  Chapters  II  and  IV. 

Let  O  be  a  convex  subset  of  R^,  and  let  {X(9),  6  6  0}  be  a  family  of 
R+ -valued  RVs. 

Definition  A.II.l.  (X(0),  6  €  0}  is  stochastically  convex  (resp.  convex  and 
increasing/ decreasing)  on  0  if  the  R^  — ►  R  mapping  6  i-»  E  f(X(6))  is  convex 
(resp.  convex  and  increasing/decreasing)  for  every  increasing  function  /  :  RhR. 
This  is  denoted  by  {X(0),  9  G  0}  €  SCX(st)  (resp.  SICX(st)/ SDCX(st)). 

When  9  is  N- valued  the  same  definition  applies  with  convexity  replaced  by 
integer  convexity.  It  is  an  easy  exercise  to  show  the  following  result  [SS.b]. 

Lemma  A.II.l.  {X(0),  9  €  0}  €  SCX(st)  (resp.  SICX(st)/SDCX(st ))  if  and 
only  if  the  complementary  distribution  function  9  1-4  P(X(9)  >  x )  is  convex  (resp. 
increasing / decreasing  convex)  for  every  x  >  0. 

First,  the  complementary  waiting  and  response  time  distributions  in  the 
M/PH/1  queue  are  shown  to  be  monotone  increasing  and  convex  in  the  system 
utilization.  Therefore,  the  waiting  and  response  times  in  the  M/PH/1  queue  are 
both  stochastically  increasing  and  convex  in  system  utilization  by  Lemma  A.II.l. 
The  results  are  then  extended  to  the  M/GI/1  queue  using  the  fact  that  the  class 
of  PH-distributions  is  dense  in  the  space  of  probability  distributions  on  [0,  00) 
[Neu.c]. 

In  this  appendix,  I  denotes  the  identity  matrix  with  appropriate  dimensions, 
while  the  m  x  m  and  the  1  x  m  row  vector  with  zero  entries  are  denoted  by  Om 
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and  Om,  respectively.  Finally,  the  notation  denotes  the  real  part  of  a  complex 
number  s. 

Waiting  and  response  time  distributions  for  the  M/PH/1  queue 

Consider  an  M/PH/1  queue  under  the  FCFS  discipline  where  the  Poisson 
arrival  process  has  rate  A  and  the  PH-type  service  time  distribution  of  order  m 
has  an  irreducible  representation  (a,  A)  with  mean  The  matrix  A  is  invertible 
so  that  absorption  into  the  state  m  + 1  from  any  initial  state  is  certain  [Neu.b,  p. 
45] .  The  corresponding  1  x  m  column  vector  of  absorption  probabilities  is  denoted 
by  A0,  i.e.,  Ae  =  —A0.  It  is  assumed  that  om+i  =  0  and  that  the  queueing  system 
is  stable,  i.e.,  p  :=  A  <  1. 

r 

Let  W(-)  be  the  stationary  waiting  time  distribution.  The  following  well- 
known  result  was  given  in  [Neu.a,  p.  181]. 

Theorem  A.II.2.  The  waiting  time  distribution  W{-)  is  PH-type  and  has  a  rep¬ 
resentation  (7,  L )  of  order  m  with 

7  =  /9  7T,  7m+l  =  1  -  P,  (A. 2. la) 

L  =  A  +  pA°  7T,  L°=(l-p)A° ,  (A. 2. lb) 

where  the  probability  vector  tt  is  uniquely  determined  by  the  relations 

7 r(A  +  A0  a)  =  0m  and  ire  =  1. 


In  an  M/GI/1  queue  in  statistical  equilibrium,  the  response  time  T  is  the  sum 
of  two  independent  RVs,  the  waiting  time  W  and  the  service  time  B.  Therefore,  in 
view  of  Theorem  A.II.2  and  of  closure  properties  of  the  PH- distributions  [Neu.b, 
p.  51],  it  is  an  easy  exercise  to  see  that  the  stationary  response  time  distribution 
T(-)  is  also  PH-type  of  order  2m  with  representation  ($,Z),  where 

(  =  (p*,(l  -  p)a)  and  Z  =  (<£_  £^“).  (A. 2.2) 
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The  following  theorem  shows  that  a  lower  order  PH-representation  with  only  m 
phases  is  in  fact  available  for  T(-). 


Theorem  A.II.3.  The  response  time  distribution  T(-)  is  PH-type  with  represen¬ 
tation  (a,L). 

Proof.  The  Laplace-Stieltjes  transform  r(s)  of  the  PH-distribution  (£,  Z)  is  given 

by 


r(s)  =  t(sl  -  Z)-1  Z°  ,  &(s)>0. 

By  making  use  of  equations  (A. 2.1)  and  (A. 2. 2),  direct  calculations  now  yield 

'( sI-L r1  (sI~L)-lL°a(sI-A)-l\  /0£ 

Om  (si -A)-'  )  \A° 


r(s)  =  (pir,(l-p)a) 


=  pn(sl  -  L)~x  L°  a(sl  -  A)-1  A0  +  (1  -  p)a(sl  -  A)-1  A0 
=  pa(sl  -  A)~1A°tt(sI  -  L)~lL°  +  a(sl  -  A)~lL° 

=  a  [p(s/  -  A)'1  A°7r  +  (si  -  A)-\sI  -  L) )  (si  -  L)~x L° 

=  a(sl  -  L)-1  L°  ,  $(s)>  0, 

thus  completing  the  proof. 

□ 

Remark  A.II.l.  The  simple  M/PH/l  queue  provides  a  building  block  in  the 
approximate  decomposition/ aggregation  algorithms  for  analyzing  queueing  models 
of  real  life  systems  [Giin.b].  In  such  an  iterative  algorithmic  analysis,  low  order 
representations  for  various  distributions  in  the  network  are  often  desirable  from  a 
computational  standpoint.  Theorem  A.II.3  serves  this  purpose  and  provides  a  PH- 
representation  for  T  with  only  half  the  dimension  of  (£,  Z).  Furthermore,  Theorem 
A.II.3  shows  that  the  representations  of  T  and  W  differ  only  in  the  way  they  are 
initialized. 

Strong  convexity  results  for  the  M/GI/1  queue 
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Theorems  A. II. 2  and  A.II.3  are  used  to  study  the  variation  of  the  RVs  W 
and  T  with  the  system  utilization  for  the  M/PH/1  queue.  The  results  are  then 
extended  to  the  M/GI/1  queue. 

The  notation  W(p)  and  T(p)  is  now  used  to  represent  the  RVs  W  and  T, 
respectively,  in  order  to  indicate  the  dependence  of  their  distributions  on  p  explic¬ 
itly.  In  particular,  monotonicity  and  convexity  of  the  functions  p  W(p,  x)  := 
P{W(p)  >  x}  and  p  T(p,x )  :=  P{T(p)  >  x}  are  established  in  Theorem  A. II. 4 
for  all  x  >  0. 

Theorem  A.II.4.  All  the  partial  derivatives  of  W(p,x)  and  T(p,x)  with  respect 
to  p  exists  and  are  positive  for  all  x  >  0.  In  particular,  the  mappings  p  t— *  W(p,x ) 
and  p  i — ►  T(p,  x)  are  both  monotone  increasing  and  convex  for  all  x  >  0. 

Equivalently,  {W(p),  p  E  [0,1)}  and  {T(p),  p  E  [0,1)}  are  both  SICX{st). 

Proof.  By  Theorems  A.II.2  and  A.II.3,  the  complementary  distribution  functions 
W (p,  x)  and  T(p,x )  are  given,  respectively,  by 

W(p,x)  =  p^A+PA°^xe  and  T{p,x)  =  ae^A+pA°n)x  e.  (A.  2.3) 

It  is  a  simple  exercise  to  see  that  if  the  matrices  dne^A+pA°n^x /dpn  have  positive 
components,  then  the  partial  derivatives  dnW(p,x)/dpn  and  dnT(p,x)/dpn  are 
all  positive  for  n  >  1. 

Since  the  matrices  A  and  A0 n  do  not  commute,  the  matrices  dne^AApA°n^  /dpn, 
n  >  1,  do  not  have  a  simple  closed  form  expression.  Therefore,  consider  the  matrix 

E(p)  :=  e0 (A+cI)+pxA0n)  _  ecX  e(A+PA°*)x  ^  x  >  q 

where  c  :=  max{— An  :  1  <  i  <  m}.  Since  both  F  :=  x(A  +  cl)  and  G  :=  xA°n 
are  positive  matrices, 

E(p)  =  eF+P0  =  f:il±lSl 

k=0 
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is  a  polynomial  in  p  with  positive  coefficient  matrices.  The  partial  derivatives 
dnE(p)/dpn  and  dne^A+pA  /dpn  are  therefore  positive  for  all  n  >  1,  and  the 
first  part  of  the  theorem  thus  follows.  The  monotonicity  and  convexity  of  the 
mappings  p  h->  W(p,  x)  and  p  >->  T(p,x )  are  now  immediate  since  dW(p,x)/dp , 
d2W(p,  x)/dp2,  dT(p,x)/dp  and  d2T(p,x)/dp2  are  all  positive  for  x  >  0. 

□ 

The  structural  result  of  Theorem  A. II. 4  is  now  extended  to  the  M/GI/1 
queue. 

Corollary  A.II.5.  IfW(p)  and  T(p)  are  the  waiting  and  response  times  in  an 
M/GI/1  queue  with  utilization  p,  then  {W(p),  p  £  [0,1)}  and  {T(p),  p  £  [0,1)} 
are  both  SICX(st). 

Proof.  The  result  follows  from  Proposition  8.2.5a  of  [Sto,  pp.  169-170]  using  the 
fact  that  the  PH-distributions  are  dense  in  the  space  of  probability  distributions 
on  [0,oo). 

□ 

Remark  A.II.2.  The  more  general  GI/GI/l  queue  is  considered  in  [SS.a]  using 
a  sample  path  approach  and  sufficient  conditions  for  sample  path  convexity  of  the 
waiting  time  in  the  parameter(s)  of  the  service  and  interarrival  times  are  given 
(see  also  [ShY]).  However,  it  does  not  seem  possible  to  establish  the  stochastic 
convexity  of  the  waiting  and  response  times  in  the  arrival  rate  from  the  results 
in  [SS.a]  or  [ShY].  Moreover,  the  stochastic  convexity  considered  here  is  stronger 
than  the  sample  path  convexity  defined  in  [SS.a]  (see  [SS.b]).  The  sample  path 
approach  allowed  the  consideration  of  a  very  general  class  of  problems  in  [FeG]. 
However,  when  specialized  to  the  M/GI/1  queue,  only  the  convexity  of  the  average 
waiting  time  in  the  arrival  rate  was  established  in  [FeG]. 

Remark  A.II.3.  In  [SS.a],  it  is  shown  for  the  GI/GI/l  queue  that,  if  the  nth 
service  time  Sn(p )  =  Xn/p  where  {Xn,  n  =  0, 1, . . .}  is  a  sequence  of  non-negative 
i.i.d  RVs,  then  the  waiting  time  of  each  job  is  stochastically  decreasing  and  convex 
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in  fi  in  the  sample  path  sense.  This  result  generalizes  a  result  established  in  [Web] 
again  by  sample  path  arguments,  namely  that  the  average  waiting  time  for  each 
job  is  convex  in  the  service  rate.  For  the  M/GI/1  queue,  the  stationary  version  of 
the  convexity  result  of  [SS.a]  in  the  stronger  SDCX(st)  sense  easily  follows  from 
Corollary  A. II. 5. 
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APPENDIX  III 


This  appendix  states  two  technical  lemmas.  Lemma  A.III.l.  is  used  in  the 
proof  of  Theorem  II.4.6,  while  Lemma  A.III.2  is  used  in  proving  Theorems  V.3.1 
and  V.5.2. 

Lemma  A.III.l.  Let  I  be  a  convex  subset  o/IR  with  jf  in  I,  and  let  f  :  1 1-»  IR 
be  a  concave  function.  If  F  :  IK  i->  IR  is  defined  by  F(x)  =  f(x k),  then 

K  1  1 

ma  x{F(x):x£lK  and  ^  xk  =  1}  =  F(— , . . . ,  — )  . 


Proof.  The  result  follows  from  the  following  simple  argument 


K  1 

F(*)  =  *£-/(*,) 

k=l  ^ 

k-1 

=  Kf(±) 

=  < . 


□ 

Lemma  A.III.2.  Let  Cti  C  IR”*',  i  =  1,2,  be  two  convex  sets.  If  the  mapping 
:  Oi  x  &2  m ►  IR  is  jointly  convex,  then  the  mapping  ^  i — >  IR  defined  by 

$(x)  =  min  ty(x,y)  ,  x  in 
y€n2 

is  also  convex. 

Proof.  The  convexity  result  follows  directly  from  Theorem  5.7  of  [Roc,  p.  38]. 

□ 
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APPENDIX  IV 


A  Proof  of  (3.3.1) 

Consider  a  stable  M/M/1  queue  with  utilization  p ,  and  let  {Tn,  n  =  1,2,...} 
be  the  Poisson  sampling  process  with  rate  v  independent  of  the  arrival  and  the 
service  processes  defining  the  M/M/1  queue.  Let  Nr  be  the  number  of  samples  in 
the  interval  [0,  r],  and  let  rjT  be  the  times  the  server  is  sampled  idle  in  this  interval. 
If  Q(t )  is  the  number  of  customers  in  the  system  at  time  t,  then 

OO 

Nr  =  £  1  (Tn  <  T) 

n=  1 

and 

Nr 

r)r  =  Y,  W(r„)  =  0) , 

n—1 

where  l(-)  is  the  indicator  function. 

Note  that  {Q(Tn),  n  =  0, 1, . . .}  with  To  =  0  is  a  Markov  chain  with  invariant 
probability  mass  7r,  i.e.,  7 r*  =  P[Q(Tn)  =  k],  k  =  0,1,...,  and  the  a.s. 

relation 

lim  -  Y]  1  (Q(Tj)  =  0)  =  7r0 

n— ►oo  n 

1=1 

thus  holds  by  the  Ergodic  theorem.  Therefore, 


Tjr 

lim  — -  =  7To  a.s. 

T— >00  Nr 


since  limr_oo  Nr  =  oo.  The  relation 


lim  =  l  —  p  a.s. 
r-» oo  Nr 


(A.4.1) 


(A. 4.2) 


now  follows  from  (A. 4.1)  and  the  PASTA  property  [Kle.a]. 

In  order  to  show  the  second  equation  in  (3.3.1),  divide  the  interval  [0,  t]  into 
n  equal  subintervals  of  length  t,  i.e.,  nt  =  r.  Let  rjk  be  the  number  of  times  the 
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server  is  sampled  idle  in  the  interval  [( k  -  l)t,  kt ),  k  =  1, . . . ,  n,  and  let  iV*  be  the 
number  of  measurements  in  this  interval.  Note  that  {Nk,  k  =  1,2,...}  is  i.i.d. 
with  E[Nk]  =  vt,  and 


n 

Vr  =Vnt  =  Yj  Vk 
=  1 


n 

and  Nr  =  Nnt  =  ^  Nk  . 

fc=i 


Therefore, 


lim 

n— ►  oo 


Y)nt 

N7t 


rjnt  n 

lim  — aTT 

n^°°  T,k=iNk  fn 


=  1  ~P 


a.s. 


by  (A. 4. 2),  so  that  the  result 


lim  -^7  =  lim  —  =  1  —  p  (A.4.3) 

n— ♦'OO  TlUt  T—+OQ  TU 


holds  since 


by  the  Law  of  Large  Numbers. 
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